Optimization of pressure swing adsorption via a trust-region filter algorithm and equilibrium theory

Optimization of pressure swing adsorption (PSA) remains a challenging task, as these are periodic dynamic systems governed by nonlinear PDE systems. This study develops an optimization strategy that incorporates reduced PSA models and also samples information from a high-fidelity ’truth’ model. The reduced model is based on equilibrium theory and can be applied to optimize a large variety of cyclic adsorption processes. The optimization is performed via a trust-region filter (TRF) method. The TRF method uses reduced models to minimize the information required from the high-fidelity adsorption model during mathematical optimization. Moreover, this method guarantees convergence to the optimum of the high-fidelity model. This approach is applied to optimize a 4 column pressure swing adsorption process, modeled by partial differential algebraic equations, for the separation of a CO 2 /CH 4 mixture. The results show that the reduced model significantly reduces the computational time of the method’s trust-region step compared to previous studies. © 2021 Published by Elsevier Ltd.


Introduction
Pressure swing adsorption (PSA) is a process commonly applied to many gas separation tasks, such as oxygen separation from air Jiang et al. (2003) , flue gas separation Xu et al. (2019) , and biogas upgrading Ferella et al. (2017) . Because of the technology's broad application in industry many studies exist on the simulation and optimization of PSA processes to maximize operation efficiencies Ding et al. (2018) and adsorbent materials. For example Dowling et al. (2012a) used a superstructure approach for simultaneous optimization of design and cycle operation, later extended by Wang et al. (2015) . An overview of the most recent advances on PSA optimization can be found in Biegler et al. (2004) Despite the corresponding improvement in mathematical modeling of the process, optimization of the PSA process is still a challenging task. Typical mathematical models for PSA processes describe the sharp adsorption fronts, cyclic behavior, and internal recycle streams by systems of stiff partial differential equations (PDEs). Spatial discretization of these models results in systems of differential algebraic equations (DAE), which are time consuming or numerically difficult to solve. The problem formulations are often ill-conditioned, and internal recycle streams and the cyclic behavior result in dense constraint Jacobians of the nonlinear programs (NLPs).
A number of approaches to PSA optimization have been demonstrated over the past three decades. One key factor of the optimization of PSA processes is the determination of the cyclic steady state (CSS), i.e., when the initial and final states of the system coincide. Many publications utilize successive substitution (i.e. Picard iteration) Ding et al. (2018) . This method iterates by integrating over the model equations, starting with the final state of the previous iteration as the new initial value. This process is repeated until the difference between the initial and final states is smaller than a predefined threshold. This method has the advantage of being numerically stable, but it exhibits slow convergence behavior. Another approach to finding a solution to a PSA model with respect to the CSS condition, is the addition of algebraic constraints. This approach has been applied to single bed PSA units by Jiang et al. (2003) and was later extended to multiple bed PSA  Jiang et al. (2004) . More recently, Tsay et al. (2018) also applied this approach to periodic moving bed chromatography. The convergence of this method is much faster, and it converges globally Biegler et al. (2004) . Simultaneously solving the fully discretized PDEs and the CSS condition has been successfully applied to set-ups with two columns. To give a few examples: Smith IV and Westerberg (1991) proposed an MINLP approach with a fully discretized PDE to optimize the operation cycle and the number of beds. The SQP method was successfully applied to optimize fully discretized PSA models by Biegler et al. (2004) and Tao et al. (2019) . Vetukuri et al. (2010) proposed a quasi-Newton method, which approximates the dense constraint Jacobian of the NLP, and could show a significant reduction in com-putational costs. Dowling et al. (2012b) formulated the PSA as an optimal control problem to determine optimal cycle times. This approach was further improved by Wang et al. (2015) , who used flux limiter schemes for reducing computational costs. Agarwal et al. (2009) developed a reduced model via proper orthogonal decomposition (POD), which was used by a trustregion filter (TRF) method for PSA optimization Agarwal and Biegler (2013) . The TRF method uses a reduced model to assist the optimization of the PDE constrained problem by reducing the evaluations of the PDE model. In Agarwal and Biegler (2013) it was noted that the computational effort of the TRF method can still be improved by choosing a tailored reduced model in the trust-region step.
Another challenge arises with the complexity of the PSA cycle. In the literature many optimization examples for 2-column PSA cycles can be found, as well as simulation studies for more complex cycles. To the authors knowledge, only Jiang et al. (2004) optimized an adsorption cycle with more than two adsorption columns. We therefore discuss the optimization of a more complex 4-column PSA cycle with a 9-step operating cycle, including 4 internal recycle streams, as opposed to the two column set-ups common in the literature.
Here we apply the TRF method as proposed in Agarwal and Biegler (2013) . We implement a different reduced model based on equilibrium theory Knaebel and Hill (1985) for the application of the TRF method to minimize the computational time of the trustregion step. Despite the increased complexity of the PSA cycle, the equilibrium model has fewer variables than the POD approach in Agarwal and Biegler (2013) .
The PSA cycle is depicted in Fig. 1 , including adsorption (AS), pressure equalization with decreasing pressure (DEQ), pressure equalization with increasing pressure (PEQ), depressurization (DP), desorption (DE) and pressurization (PR). The cycle includes four recycle material streams: two pressure equalization streams, one purge stream, and one stream to repressurize the column prior to the adsorption step. The pressure is constant during AS and DE. The product gas methane is produced at high pressure during the adsorption step and is partially returned to the column during PR. The off-gas is produced at low pressure during DP II and DE. Gas is exchanged during the pressure equalization steps, from DEQ I to PEQ I and from DEQ II to PEQ II. Furthermore, gas from the depressurization DP I is used to purge the column at DE and a fraction of the AS product gas repressurizes the column at PR. We determine the final product gas concentration at cyclic steady state (CSS). As an example, we model the separation of methane and carbon dioxide, as commonly applied for biogas upgrading. The goal is the purification of the light component methane from the mixture via adsorbents such as carbon molecular sieves with a focus on product purity, recovery and process efficiency. As a second example, we consider the separation of a ternary mixture, including hydrogen. This application may arise with the production of synthetic natural gas from hydrogen and carbon dioxide Uebbing et al. (2020) . As hydrogen is almost inert on the adsorbent, hydrogen remains in the product methane. Only carbon dioxide is separated from the mixture as an impurity.
We note that these applications are examples to model the separation process. The method and the proposed reduced model can be applied to any other adsorption separation task by applying different model parameters and boundary conditions. The next section presents the detailed PDAE model for the PSA system along with the reduced model based on equilibrium assumptions. Section 3 introduces the trust-region filter algorithm and its adaptation to high-fidelity optimization with reduced models. Results are presented in Section 4 for two optimization case studies for the PSA system. Section 5 concludes the paper with a summary and directions for future work.

Modeling
In this section we introduce the high-fidelity and reduced PSA models: In Section 2.1 we describe the PSA model, which we aim to optimize in this work. The high-fidelity PSA model is described by a set of partial differential algebraic equations (PDAEs) with changing boundary conditions, according to the cycle configuration in Fig. 1 . Optimization of this PDAE model is computationally and numerically quite challenging. Much simpler representations for PSA models have been derived in the past. One example is the equilibrium model by Knaebel et al. Knaebel and Hill (1985) , which we introduce in Section 2.2 . The equilibrium theory allows for reformulation of the PDAEs to algebraic equations, which are much easier to evaluate and optimize.

PSA modeling via PDAEs
The PSA is described by a system of partial differential and algebraic equations assuming ideal gas behavior, no axial pressure gradient, no accumulation in the shockwaves, non-isothermal adsorption behavior, adsorption according to the linear driving force model (LDF), as used by Park et al. (20 0 0) , and Langmuir-type equilibrium isotherms Canevesi et al. (2018) . The mass balances do not include diffusion terms to avoid smearing of the steep adsorption fronts. We denote u the interstitial velocity, p pressure, T temperature, T w wall temperature, q i the amount adsorbed of component i ∈ { CH 4 , CO 2 , H 2 } , and q * i the amount adsorbed at adsorption equilibrium. Furthermore, let us define f s = 1 − ρ s , where denotes the bed void fraction and ρ s the density of the adsorbent.
A detailed list of the notation of the PDAE model can be found at the end of this work. The partial mass balances of the components i ∈ { CH 4 , CO 2 } are given by The sum of mole fractions determines the partial pressure of hydrogen. The energy balances (3) and (4) model the heat transfer between the adsorbent and the gas, as well as the gas and the column wall.
(c g C + c s f s ) ∂T ∂t The adsorption equilibrium for component i ∈ { CH 4 , CO 2 } is given as a multi-component extension of the multi-side Langmuir isotherm.
Parameters for the adsorption of methane and carbon dioxide of the adsorbent CMS-KP 407, are taken from Canevesi et al. (2018) .
Hydrogen is assumed to behave as a non-adsorbing gas component. The pressure in the column changes according to Here α ≥ 0 determines the speed at which the pressure of the column changes. A large value for α denotes fast pressurization or depressurization, while a small α denotes slow pressurization or depressurization. The total mass balance ∂C ∂t in combination with the ideal gas law C = p/ (RT ) gives We use (3), (7) , and (8) to eliminate the derivatives in time and get to determine the interstitial velocity u, where The PDAE system is discretized in space via a finite volume method to get a system of differential algebraic equations (DAEs). We consider the cyclic operation with 4 columns and 9 steps Uebbing et al. (2019) , shown in Fig. 1 , including pressure equalization steps between the columns.
During the different steps of the adsorption cycle different boundary conditions apply. The boundary conditions show the connections of the four internal recycle streams. The gas concentrations during the desorption step y i,purge are calculated from the average product gas of DP I. The interstitial velocity u purge is calculated such that the complete product of the step DP I is fed back to the column during the DE step. In the same way, the feed concentrations and velocities of the pressurizing PEQ I and PEQ II steps are connected to the respective depressurizing steps. Here, the parameter α is free to avoid over-determination of the system variables. For the other PSA cycle steps the value of α is fixed; to α = 0 . 01 during the DP I step, to α = 0 . 5 for the DP II step and α = 0 . 1 otherwise. This results in the rate of pressure change adapting to the gas flow into the column. The composition of the gas fed to the column during the pressurization step PR is equal to the product gas of the adsorption step AS, i.e. a part of the product gas is used to re-pressurize the column. For the pressure parameters p AS > p PEQI > p DPI > p PEQII > p DE holds. The intermediate pressures are determined as linear functions of the ad-and desorption pressure. The time of the respective PSA step is determined as t s · t f , where t f denotes the total cycle time. Table 1 shows the boundary conditions of the DAE system at the time of the PSA step.
The pressure changes according to p end (t) which is constant during each step of the PSA cycle. The value assigned to p end in each PSA step is shown in Table 1 . Fig. 2 shows an example of the pressure in the adsorption column over time.
Semi-discretization of the PDAE model results in a DAE system with differential states termined by the discretized PDEs (1), (7), (3), (4) , and (8) . The al- , u , α PEQI , α PEQII ] are given by (2), (5), (9) , and the additional boundary conditions in Table 1 . Note here that x PSA and y PSA do not correspond to the variables x and y of the TRF method introduced in Section 3 . The system is difficult to solve directly via full discretization to an NLP. The sharp adsorption fronts result in stiff differential equations and internal recycle streams during the pressure equalization steps complicate the solving process further. Therefore, we calculate the CSS via the stable method of successive substitution. The Table 1 PSA boundary conditions and pressure change The boundary conditions determine the velocity, concentration, and temperature of gas flowing into the column. If u = 0 holds at the end of a column, Neumann boundary conditions of the form ∂y ∂t = 0 hold for the mole fractions y i and gas temperature T . The pressure in the column changes towards the value of p end according to (8) .

PSA step
3/20 1/20 successive substitution terminates when the CSS error is below a predefined threshold denotes the states of the PSA model at t = 0 . Implementation of a direct method to reach the CSS, as mentioned in Section 1 , can further improve the proposed method.

PSA modeling via equilibrium theory
The PSA model introduced by Knaebel and Hill (1985) is able to represent the dynamic behavior within the PSA columns without the need for partial differential equations. Additional assumptions on the PSA, namely isothermal behavior, amount of gas adsorbed is always at adsorption equilibrium, linear adsorption isotherms and a binary gas mixture, facilitate the PDAE model. The result is a model, for which it is possible to find analytical solutions. We give a brief overview over the model equations and refer to Knaebel and Hill (1985) for the detailed derivation.
The different steps of the PSA operation cycle shown in Fig. 1 can be classified into two types: PSA steps with changing pressure (DEQ I, DP I, DEQ II, DP II, PEQ II, PEQ I, PR) and PSA steps with constant pressure (AS, DE). Let denote the linear adsorption isotherm, where q A is the adsorbed amount of component A . Furthermore, let A and B (i.e., CO 2 and CH 4 ) be the heavy and light components, respectively, z ∈ [0 , 1] is the dimensionless position in the adsorption bed, ≤ 1 . We also consider the partial mass balance in terms of partial pressure of A and the total mass balance Using the assumption of linear adsorption, we reformulate these mass balances to get Applying the method of characteristics to Eq. (10) gives The model includes the position of shockwaves and waves within the adsorption column. Shockwaves appear, if a gas mixture moves towards a gas mixture of higher concentration of light component. If a gas mixture moves towards a gas mixture of lower concentration of light component, a wave occurs. A shockwave occurs during the steps AS, DEQ I, DP I, and DEQ II, a wave appears during the steps DE, PEQ II, PEQ I, and PR. In these cases, the reference points, to determine the state of the column, are chosen from the same side of the shockwave or wave. Under the assumption that there is no accumulation at the shockwave, holds for the interstitial velocity of the shockwave u s . The subscript L denotes a value directly in front (lead) of the shockwave and T denotes a value directly after (trail). If the pressure in the adsorption bed changes, u = 0 holds at one of the ends of the adsorption bed. We assume without loss of generality that u = 0 at z = 0 . Then integration of (11) gives and from (13) and (15) we get where the subscript 0 denotes a reference value on the same characteristic. From (14) follows By considering u and u s as functions over pressure, we get a system of algebraic equations to describe the PSA steps with changing pressure.
If the pressure is constant, i.e., ∂ p ∂t = 0 , integration of Eq. (11) results in From (13) holds for the shockwave.
To determine the initial conditions of the PSA column in the next PSA step, we consider linear approximations of the previous state, still preserving the position of shockwaves. If a shockwave is in the column during a change in PSA steps, the states above and below the shockwave are approximated separately. The CSS condition is added to the model as an algebraic constraint. Fig. 3 shows an example of the behavior of the characteristics for the different boundary conditions of the PSA cycle.
In Fig. 3 A the characteristics represent the AS step (and the shockwave), with feed gas entering the column from the top and product gas leaving the bottom at constant pressure. Fig. 3 B and C show different depressurization steps (and a shockwave) with gas leaving one end of the column. In Fig. 3 B the top of the column is closed, gas leaves at the bottom, representing the steps DEQ I, DP I, and DEQ II. Fig. 3 C shows the reversed case of step DP II, during which gas is leaving the column at the top. The DE step is represented by Fig. 3 D (along with a wave). During the DE step purge gas enters the column from the bottom and off-gas leaves the column at the top at constant pressure. Repressurization of the column is done via gas entering the column from the bottom, as is shown in Fig. 3 E. Fig. 3 E corresponds to the steps PEQ I, PEQ II, and PR.

Trust-region filter method
In this section we introduce the TRF method for optimization of nonlinear optimization problems as developed by Eason (2018) . We first present the main idea, before highlighting the additional parts of the algorithm. Finally, we summarize the assumptions needed for convergence of the TRF method. For a detailed proof of the convergence of the TRF method we refer to ( Eason, 2018; Eason and Biegler, 2019 ).

Main idea of the TRF method
The TRF method finds the optimal solution of a nonlinear problem of the form We denote constraints g and h as glass box constraints, as they are easy to evaluate and differentiate. The constraints d, on the other hand, are very time-consuming and numerically difficult to evaluate. The TRF method reduces the number of calls to d(w ) : R n w −→ R n y during the optimization, by replacing d with a local surrogate model. A local surrogate model r k (w ) : R n w −→ R n y , i.e., a reduced model (RM) at the current iterate, replaces d(w ) , which we call the truth model (TM), in each iteration of the algorithm. Instead of (16) a series of subproblems is solved within a respective trust-region k .
For the convergence of the TRF method, r k (w ) must be κ-fully linear in the trust-region at each iteration. This condition assures that for k −→ 0 the difference between the two models and their sensitivities within the trust-region converge to zero. If the model sensitivities ∇d(w k ) are known, we get a κ-fully linear reduced model r k (w ) from any sufficiently smooth model ˆ r k (w ) : R n w −→ R n y by applying the the first order correction, defined by

Additional strategies and pseudocode of the TRF algorithm
The TRF method includes additional strategies to handle infeasible subproblems and to determine conditions for termination. The algorithm's pseudocode is shown in Algorithm 1 .
To make sure that a feasible point to (17) close to the trustregion center exists, where the approximation of the κ-fully linear RM is more reliable, the TRF method solves a series of compatible problems (17) .
Definition 3.2. Compatibility: The trust-region step from (17) is If (17) is not compatible, the algorithm enters a restoration phase to create a new iterate x k +1 and a new RM r k +1 (w ) , which results in a compatible subproblem (TRSP k +1 ). For the restoration phase to be successful it is sufficient to find a feasible point to (16) .
If x k +1 is a feasible point to (16) and r k +1 (w ) was created via the correction (19) , the new subproblem (TRSP k +1 ) is guaranteed to be compatible, as the trust-region center x k +1 is a feasible point. For example, a feasible point to (16) can be found by repeatedly solving the optimization problem  1 , 2  Optimize trust-region step from TRSP k wrt. r k (w ) , k and obtain x * k . 10: if ˆ x k acceptable to the filter then 11: x k +1 ← x * k 12: Restoration: Find x k +1 , which is feasible for ( ?? ) via iteration over the tear stream w i and choosing w i +1 = w * from the optimal solution.
To get an indicator of how close the current iterate is to an optimal point, we consider the criticality measure. We now define x c,k to be a feasible point for (17) Then the criticality measure χ k is given by The criticality measure χ k goes to zero, if the iterate x k approaches a KKT-point of (17) without the trust-region constraint.
Because the RM is κ-fully linear, the error of the RM approaches zero for k −→ 0 , and a KKT-point of (17) without the trust-region constraint approaches a KKT-point of (16) . Hence, if the criticality measure is small with respect to the trust-region radius, the trust-region radius is reduced and the TRF method continues until k approaches zero. On the other hand, shrinking k to 0 is not needed if the RM is generated via the first order correction (19) , ∇r k (w k ) = ∇d(w k ) holds, and χ k = 0 indicates a KKT-point of (16) .
In this case. optimality holds, even if k is large.
The TRF method furthermore includes a filter check. A filter is defined as the set , if sufficient progress was made to improve feasibility or objective of the previous iterates. If a new step is rejected by the filter, the iteration continues with x k +1 = x k and a reduced trust-region radius. If a step is accepted, the switching condition is checked. If (20) holds, the iterate is an f-type step. In this case, the new iterate is accepted and the trust-region radius k is increased. If (20) does not hold, the iterate is a θ -type step. The previous iterate ( f k , θ k ) is added to the filter, and the trust-region radius is changed according to An exception is made if ρ k < 0 holds. In this case a new step is rejected, despite making small progress in f k .

Convergence
Eason (2018) showed that the TRF method converges to a first order KKT point of (16) , given the following assumptions hold: A1 The functions f, g, h, and d defining (16) are twicecontinuously differentiable. A2 The problem domain is closed and bounded. A3 The Mangasarian-Fromovitz constraint qualification (MFCQ) holds for (16) at all limit points of the TRF iteration.

A4
The reduced model is κ-fully linear, twice-continuous differentiable and the second derivatives are uniformly bounded. A5 The solution ˆ x k of (17) reduces the objective function value according to the fraction of Cauchy decrease for a κ t > 0 and a bounded sequence β k > 1 . In other words, the solver used for optimizing the trust-region step must make sufficient progress in relation to the criticality measure evaluated in x c,k . In practice this condition is fulfilled, by using an NLP optimization strategy to solve (17) , which is initialized in the feasible point x c,k .

A6
The condition x c,k − x k ≤ κ u θ k holds for small θ k < δ and a κ u > 0 . If r k (w k ) = d(w k ) holds, as is the case if the first order correction (19) is used to generate the RM, this condition is fulfilled.
Eason shows that the TRF method will create a subsequence { k i } with compatible trust-region steps (TRSP k i ) and lim i −→∞ where x * is a KKT-point of (16) . For the convergence proof we refer to Eason (2018) .

Trust-region radius
As shown by Yoshio and Biegler (2021) , the trust-region radius in the subproblem (17) must not necessarily extend to all model variables x . Instead it can be formulated in terms of the degrees of freedom alone. We partition are the degrees of freedom and ˆ As long as the model sensitivities are non-singular, the trust-region constraint on the degrees of freedom propagates to the remaining variables according to For a detailed proof we refer to Yoshio and Biegler (2021) . We therefore can rewrite the trust-region radius as In the current work, the degrees of freedom correspond to x = w . We apply both strategies, the full trust-region radius and the trustregion radius regarding the degrees of freedom, and compare the results.

Regularity and feasibility of the trust-region step
To ensure that MFCQ holds, which is required in the limit point of the iteration to ensure convergence according to assumption A3 , for all x ∈ R n , one can introduce artificial variables and 1 penalties, and rewrite (16) as: where e h ∈ R n h and e g ∈ R n h + n y are vectors with elements of 1 and E is a scaling matrix. Note here that the artificial variables only need to be added to constraints which may violate the MFCQ. The corresponding trust-region subproblem has a feasible solution with x q, j = max (0 , g j (x )) , x p ⊥ x n , x p,i + x n,i = | [ h (x ) , (y − r k (w )) ] i | and the compatibility check can be skipped. In this case, we compute the infeasibility measure as and enter the restoration phase if k ≤ tol and θ k > θ tol .

Application of the TRF method to PSA optimization
We show here, how we optimize the PDAE model from Section 2.1 with the TRF method and how the model based on equilibrium theory from Section 2.2 is used as a local surrogate model. In Section 3.5.1 we define the function d(w ) , which is used to apply the TRF method. The TM d(w ) represents the correlation between the columns design and cycle operation, given by the variable w, and the product gas flow and concentrations of the PDE model at cyclic steady state (CSS), denoted y . In Section 3.5.2 we show the corresponding reduced model r(w ) , which calculates the correlation between design and product via a set of algebraic equations derived from equilibrium theory. Finally, we discuss the calculation of derivatives of the TM, which are needed to apply the first order correction.

The truth model d(w )
To optimize the separation performance of the PSA via the TRF method, we need to introduce the variables w, y, and the TM function d(w ) . The degrees of freedom of the PSA model are the adsorption pressure p AS , desorption pressure p DE , column diameter R i , column length L, cycle time t f , and a fraction of product gas fed back to the column during the PR step B f . With proper scaling of the variables, we define w = [ p AS , p DE , R i , L, t f , B f ] . The model response is given by the gas flow rates of the product stream which is the flow rate of the gas leaving adsorption column during the adsorption step AS. This implies that for each function call d(w ) the CSS of the DAE system must be evaluated.

Building a reduced model
The reduced model r(w ) , which is needed for the TRF method, has the same output stream r k (w ) ] as the TM.
The reduced model is calculated based on dimensionless pressure and time. The parameters R i , L, and t f are scaling the interstitial velocity of the gas flowing in and out of the column. In addition to the inputs w for the truth model, the equilibrium model has additional degrees of freedom, which are the adsorption parameters k A , k B as well as the concentration of the gas entering the column during DE and PR, denoted y DE and y PR . We allow for different values of the adsorption parameters in various step of the PSA process, . . ] and let p RM denote the additional DOFs p RM = [ k A , k B , y DE , y PR ] . We use these additional DOFs to derive a local reduced model r k (w ) from the equilibrium model at a current iterate w k of the TRF algorithm: Let r(w, p RM ) denote the model response of the equilibrium model, which is the product gas flow rates of the separation at CSS. We minimize the error of the equilibrium model to the TM response d(w k ) at the current iterate according to and define ˆ r k (w ) = r(w, p * k ) . To guarantee convergence of the TRF Eq. (18) must hold for the reduced model at the current iterate w k . We can assure that these conditions hold, by applying the first order correction (19) .

Derivatives
To apply the first order correction (19) we need to calculate the sensitivities ∇d(w k ) of the TM at CSS. One option to calculate the sensitivities is the finite difference approach, which is simple to implement, but has several disadvantages in practice. Firstly, the evaluation of the finite differences is very time consuming. It requires multiple function evaluations of the TM, each of which include the calculation of the CSS via successive substitution. Furthermore, the successive substitution calculates the CSS only up to a predefined tolerance || An alternative is to consider d as a function of the input w and x PSA (0 , w ) at CSS, which is implicitly depending on w via the CSS where the partial derivatives ∂d(w,x PSA (0 ,w )) ∂w and ∂d(w,x PSA (0 ,w )) ∂x PSA (0 ,w ) are the backward sensitivities of the PSA model equations. Furthermore, we apply the implicit function theorem to get from the CSS equation. We then use automatic differentiation of CasADi Andersson et al. (2014) to calculate (25) and the backwards sensitivities.

Comparison of TM and RM
To get an impression of the accuracy of the RM, we compare the results at the reference point w re f for separation of a binary gas mixture. The values of the reference point are shown in Table 2 .
For the comparison, we consider the separation of a binary mixture of CO 2 and CH 4 . Fig. 4 shows the mole fraction of CH 4 in the gas phase of one column over bed length (ordinate) and time (abscissa), starting with the adsorption step AS. Notable is a shift in time of the desorption of CO 2 from the adsorbent, which is indicated by the vertical dark blue area in this figure at times 0.5 to 0.7. Also, while the TM shows a monotonic decrease of CH 4 concentration over space, the RM has increased CH 4 concentrations near the end of the column. Both of these phenomena can be explained by the adsorption kinetics of the models. The ad-and desorption happens instantaneously in the equilibrium RM, while the TM has adsorption kinetics that slow the ad-and desorption. Fig. 5 shows the divergence of amount adsorbed and the adsorption equilibrium in the column of the TM. The strong adsorption of CO 2 during the desorption step, at times 0.5 to 0.7 on the dimensionless time scale, is clearly visible, causing the aforementioned shift in time. The difference in the amount adsorbed and the amount adsorbed at equilibrium is higher for CH 4 , which is the light component, because of the faster adsorption kinetic of CO 2 .
Important for the speed of convergence of the TRF method is not the accurate representation of the states within the column, but an accurate representation of the TM model response, d(w ) , close to the reference point at which the RM was created. Fig. 6 shows the relative error of the model response with respect to the distance to the reference point w re f . This figure shows that the error is small close to the reference point, as desired, and increases linearly with greater distance.

Optimization of PSA processes via TRF method
We apply the TRF algorithm to optimize the PSA model with respect to recovery and purity of the product gas CH 4 . The trust region step subproblem was solved using IPOPT Wächter and Biegler (2006) and MA57 HSL (0000) . The variables w, y, and the function d(w ) are as defined as in Section 3.5 . We introduce additional variables, v = [ v p , v r ] , which represent the purity and recovery of the product methane and define the optimization problem The objective corresponds to finding a Pareto optimal point with respect to product purity and recovery.
The superscript f eed refers to the mole flow rate of the feed gas. We choose a CSS tolerance of to assure that the error of the model evaluation decrease with the trust-region radius. The results we show here have a trust- Fig. 5. The difference q * − q between amount adsorbed q and amount adsorbed at equilibrium q * in mmol/g. The plots show the difference of the amount adsorbed q and the amount adsorbed equilibrium q * with parameter w re f and feed y CO 2 = 0 . 4 , y CH 4 = 0 . 6 for the components CH 4 and CO 2 in mol/g.  (19) . After fitting the parameters of the equilibrium model to the model response of the TM at reference point y TM = d(w ) , the resulting RM was evaluated at different inputs w on a grid around w re f . We define y RM = r re f (w ) . The x-axis shows the difference of the input w to the reference point w re f . At w = 0 the error of the parameter fitting, which was used for building the RM, is shown. The plot shows the increasing error of the model with greater distance from the reference point w re f . region radius of 10 −4 to 10 −6 upon termination, which implies a CSS tolerance of 10 −7 to 10 −9 . The parameters of the TRF algorithm for the different cases we consider in the following are shown in Table 3 . We also validate the assumptions for convergence of the TRF method introduced in Section 3.3 . Assumption A2 is fulfilled by having box constraints for the problem variables. By using the first order correction (19) and proper formulation and scaling of the RM, Assumptions A4 and A6 hold. Initializing the optimization of the subproblem in the trust-region step with x c,k and returning a Objective is to find a Pareto optimal point with respect to purity and recovery. The trust-region radius is calculated with respect to all variables x . local optimum results in sufficient progress with respect to x c,k to fulfill A5. During the optimization run, we check, whether the linear independence constraint qualification (LICQ) holds for (16) . The LICQ holds at every iterate of our optimization runs and implies the MFCQ, which is needed for A3. Finally, Assumption A1 may not always hold, because the model response d(w ) and sensitivities ∇d(w ) include the error of the CSS calculation and the integration of the discretized model DAEs. During the optimization this is noticeable by oscillations in the objective value, feasibility measure and criticality measure. The error of the model response d(w ) was empirically determined to be around 10 · CSS tol , thus for a CSS tolerance of 10 −8 we have a corresponding model error of around 10 −7 . This error propagates to the sensitivities, because for the application of the implicit function theorem (25) we assume to be at the CSS. Because the oscillations occur at a larger order of magnitude however, the errors are mostly caused by the integrator. We terminate the TRF method here, when θ k , χ k ≤ 10 −5 .
In the following, we separate a 4:6 binary mixture of CO 2 and CH 4 . Fig. 7 shows the development of objective f k = f (x k ) , infea- sibility measure θ k = y k − d(w k ) and criticality measure χ k over time, where the trust-region radius is calculated with respect to all variables. Termination occurs after 72 iterations. We validated the calculation of derivatives described in Section 3.6 with a straightforward Finite Differences approach. Qualitatively, the calculated derivatives were identical. However, as expected the unavoidable numerical noise in the Finite Differences led to an increase in the number of iterations (roughly 300). The number of iterations needed and the total time spend in the calculations of the trust-region step are summarized in Table 4 as case A (TR: full). In the previous work of Agarwal and Biegler (2013) a reduced model based on proper orthogonal decomposition (POD) was used for the optimization of a 2 column PSA cycle, which has 52247 algebraic variables. For a direct comparison we note that the study of Agarwal differs from ours in multiple aspects, such as the PSA set-up, the number of modeled columns, the cycle configuration, and the direct determination of the CSS via Newton method. Considering these differences, we can make the following observations; our approach requires more iterations, presumably due to the lower model accuracy of the RM and the more complex PSA set-up. However, the total accumulated time spend in the trust-region step is significantly lower. We attribute the reduction in computational time to the smaller size of the RM. The proposed RM based on equilibrium theory has a total of 67 variables for our 4-column, 9-step configuration and required a total of 5.5 seconds to solve over 72 iterations. If the trust-region radius is only calculated with respect to the degrees of freedom, termination occurs after 55 iterations. The results are shown in Fig. 8 and Table 4 as case A (TR: DoFs).
As a second case, we consider the calculation of the trust-region radius with respect to the degrees of freedom according to (21) . The optimization results are shown in Fig. 8 and in Table 4 as case A (TR: DOFs). Table 4 shows that fewer iterations are needed in this case.
The equilibrium model which we apply as the reduced model has the drawback of only modeling binary mixtures in the gas phase. To analyze if the optimization of a ternary mixture is possible nonetheless, we optimize the separation of a 1:1:1 mixture of CO 2 , CH 4 and H 2 , which we denote as case B (H 2 , TR: DoFs) in Table 4 . The hydrogen in the gas phase is here approximated only by the first order correction term. In this case the criticality measure did not reach the threshold of χ tol = 10 −5 with the adaptive CSS tolerance (26) even after 200 iterations. We attribute this to the lower accuracy of the reduced model in this case, and chose a fixed CSS tolerance of CSS tol = 10 −8 instead of using (26) . While this increases the overall computational costs, it converges successfully and improves the reliability of the function evaluations and corresponding gradients. The result is shown in Figure 9 . Objective f k , infeasibility measure θ k , and criticality measure χ k over the number of iterations (trust-region radius: DoFs). As in Fig. 7, but now the trust-region radius is calculated with respect to the degrees of freedom. Fig. 9. Objective f k , infeasibility measure θ k , and criticality measure χ k over the number of iterations (mixture with hydrogen, trust-region radius: DoFs). As in Fig.  8, but a ternary mixture with hydrogen is the feed gas to the PSA columns.

Optimizing PSA work demand
In the following we use the introduced algorithm and the reduced model to the application of biogas upgrading. Carbon dioxide must be removed from the product gas of anaerobic digestion, which is a gas mixture of CO 2 and CH 4 , prior to feed to the gas grid. Typically, product gas from anaerobic digestion contains 50 − 70 mol.-% of CH 4 Mohseni et al. (2012) . We assume here that the mixture contains 60 mol.-% CH 4 and 40 mol.-% CO 2 . The purity of the product gas must meet specifications of the local gas grid, i.e., 95 mol.-% of CH 4 after the separation, according to Adler et al. (2014) .
We formulate a new objective which we aim to minimize, where is the specific reversible isothermal work demand for the gas compression at T amb = 298 . 15 K, η comp = 0 . 8 is the efficiency of the compressor, R = 8 . 3145 · 10 −3 kJ/mol/K is the gas constant, and e CH 4 = 831 . 9 kJ/mol is the chemical exergy of the product methane. The objective reflects the chemical exergy of the produced methane reduced by the work demand for the purification in kJ/mol. Furthermore, we add the constraint v p ≥ 0 . 95 (27) to the model to ensure the desired product quality in terms of mole fraction. We use this optimization problem to compare the two approaches to ensure feasibility of the trust-region step, which we introduced in Section 3.4 : On the one hand, we use the compatibility check to test the feasibility of the trust-region step. If the trust-region step is infeasible, the restoration phase is called. We refer to this approach as case C (Compatibility Check) in the results of Table 4 . Alternatively, the trust-region step is always feasible if artificial variables are added as described in Section 3.4 . We call this approach case C (Artificial Variables). A fixed CSS tolerance of CSS tol = 10 −8 was applied as for case B. The parameters of the TRF algorithm are shown in Table 3 . Fig. 10 shows the optimization of the PSA model with the compatibility check until termination with θ k = χ k = 10 −5 . The computational time is listed in Table 4 . The iteration terminates after 94 iterations. The final optimization leads to a work demand of 12.1 kJ/mol for the purification of raw biogas. According to Bauer et al. (2013) a work demand of 0.15 to 0.3 kWh/Nm 3 is typical in the industrial application of PSA processes for biogas upgrading. Assuming a mole density of 4 4.4 4 mol/Nm 3 for the raw biogas, this corresponds to 12.2 to 24.3 kJ/mol. The result of our optimization is therefore slightly below the minimum level of technology reported in the industry. With the sensitivities of the true model we can determine the KKT multiplier of the purity constraint (27) at the optimal solution. This multiplier has a value of 3.39, which indicates that relaxing the product purity to 94 mol.-% decreases the objective value by approximately 33.9 kJ/mol. This is mostly caused by an increase of product recovery v r .
Alternatively, Fig. 11 shows the optimization of the PSA model with artificial variables. In this case we stop the optimization at iteration 73, when a criticality measure of χ k = 1 . 02 · 10 −5 is reached. To reach a point with χ k ≤ 10 −5 a total of 209 iterations are needed here, despite that almost no progress is made after iteration 73 (relative difference of the objectives at iteration 73 and 209 below 0.0 0 01 %). Here, the criticality measure fluctuates strongly without noticeable progress on either objective or feasibility. Initially the error in the CSS condition, i.e. a too large CSS tolerance, was suspected to cause the very slow convergence in Fig. 10. Objective f k , infeasibility measure θ k , and criticality measure χ k over the number of iterations (with Compatibility Check). As in Fig. 7, but the objective is to maximize the units efficiency with respect to a product purity of at least 95 %. The feasibility of the trust-region step is enforced by performing the compatibility check and calling the restoration phase if necessary. Fig. 11. Objective f k , infeasibility measure θ k , and criticality measure χ k over the number of iterations (with artificial variables). The difference to the case shown in Fig. 10 is that the trust-region step is always feasible because of artificial variables, which are penalized in the objective. The restoration phase is called if an iteration is not feasible or optimal, but the trust-region size is below the predefined threshold. this case. However, no change in this behavior was observed when evaluating the true model with the very tight CSS tolerance of 10 −8 used here instead of the adaptive CSS tolerance (26) . Thus, errors of the integrator cause the iteration to stall, not the CSS tolerance. The final objective value is shown in Table 4 as case C (Artificial Variables). The value is slightly larger than the final objective of the comparative case C (Compatibility Check). We attribute this to the fact that we apply a local solver to a non-convex NLP, which can result in termination at a different local optimum.

Conclusions
We propose a reduced model based on equilibrium theory, suitable for optimization of PSA processes within the trust-region fil-ter (TRF) algorithm. Our reduced model results in a significant reduction in computational cost of the trust-region step, even for complex PSA cycles, over a comparative study by Agarwal and Biegler (2013) . We attribute the reduction in computational time to the reduced number of variables of the model. Nevertheless, by sampling the truth model (TM), the TRF method converges to the optimal solution for the truth model.
The reduced model was applied successfully to optimize a 9step 4-column PSA process. The optimization required more iterations than the comparative study, yet the overall computational time spent in the trust-region steps was reduced significantly, from 4 to 6 seconds as opposed to 4896 seconds with a POD approach.
The trust-region filter method we applied required fewer iterations until termination if the trust-region is calculated only in the degrees of freedom (55 iterations) as opposed to all variables (72 iterations).
The solution of the reduced model based on equilibrium theory differs from the PDAE (TM) solution within the column. We identified that assuming adsorption to be at equilibrium is the main reason for the differences between the two models. We expect processes with fast ad-and desorption kinetics to be represented more accurately by the reduced model which could reduce the number of iterations needed for convergence.
The bottleneck of the proposed method is the calculation of the cyclic steady state (CSS). An implementation of simultaneous methods for CSS calculations, as previously applied in the literature would shorten runtime and improve the accuracy of calls to the PSA model. See for example Jiang et al. (2004) and Tsay et al. (2018) , as introduced in Section 1 .

Declaration of Competing Interest
This research was supported by the Center of Dynamic Systems (CDS), funded by the EU-program ERDF ( European Regional Development Fund ). The funding source had no involvement in the collection, analysis and interpretation of data, in the writing of the report or in the decision to submit the article for publication.