Rebello, C. M. et al.

Optimizing CO_2 Capture in Pressure Swing Adsorption Units: A Deep Neural Network Approach with Optimality Evaluation and Operating Maps for Decision-Making

1]Carine Menezes Rebello[

1] Idelfonso B. R. Nogueira[

[1]organization=Department of Chemical Engineering,
addressline=Norwegian University of Science and Technology,

[1]

[cor1]Corresponding author

This study presents a methodology for surrogate optimization of cyclic adsorption processes, focusing on enhancing Pressure Swing Adsorption units for carbon dioxide ($CO_2$) capture. We developed and implemented a multiple-input, single-output (MISO) framework comprising two deep neural network (DNN) models, predicting key process performance indicators. These models were then integrated into an optimization framework, leveraging particle swarm optimization (PSO) and statistical analysis to generate a comprehensive Pareto front representation. This approach delineated feasible operational regions (FORs) and highlighted the spectrum of optimal decision-making scenarios. A key aspect of our methodology was the evaluation of optimization effectiveness. This was accomplished by testing decision variables derived from the Pareto front against a phenomenological model, affirming the surrogate model's reliability. Subsequently, the study delved into analyzing the feasible operational domains of these decision variables. A detailed correlation map was constructed to elucidate the interplay between these variables, thereby uncovering the most impactful factors influencing process behavior. The study offers a practical, insightful operational map that aids operators in pinpointing the optimal process location and prioritizing specific operational goals.

* Enhancing PSA efficiency with AI-based multi-objective optimization.

* The validation of optimization through artificial neural networks is important to assess the reliability of the optimized process.

* Understanding the viable operating region (FOR) helps identify the practical limits within which the system operates effectively and efficiently.

Pressure swing adsorption Multiobjective optimization Surrogate model Feasible operation region Optimally

$K$Compact subset of $\mathbb{X}$
$H$Generic prediction function or Enthalpy [J/mol]
$A$Area of column [m$^2$]
$\omega$Prediction uncertainty
$ln$Natural logarithm
$\mathbf{\theta}$Decision variables in Fisher test
$\lambda$Vector of Lagrange multipliers in Fisher test
$n_{\mathbf{\theta}}$Number of decision variables
$\alpha$Confidence level of the Fisher–Snedecor test
$\varepsilon$Predetermined tolerance
$k$Order of the function or order parameters

$L$Fisher test function or length of column [m]
$P $ Pressure [Pa]
$P_0 $ Adsorption pressure [Pa]
$P^* $ Dimensionless pressure
$T $ Temperature [K]
$T_0 $ Feed temperature [K]

$x_i $ Dimensionless molar loading of component $i$ in the solid phase
$q_i $ Molar loading of component $i$ in the solid phase [mol/kg]
$q_{s0}$ Molar loading scaling factor [mol/kg]
$\nu_0 $ Velocity scaling factor [m/s]
$\nu_z $ Superficial gas velocity [m/s]
$\nu^*_z $ Dimensionless superficial gas velocity
$Z $ Dimensionless length coordinates
$z $ Bed length coordinate [m]
$x_i^* $Dimensionless equilibrium molar loading of component
$i$ in the solid phase
$q_i*$ Equilibrium molar loading of component $i$ in the solid
phase [mol/kg]
$R $ Universal gas constant [J/(mol.K)]
$D$Inner diameter of column [m]
$D_m$Molecular diffusivity of $CO_2-N_2$ mixture [m$^2$/s]
$D_L$Axial dispersion coefficient [m$^2$/s]
$r_p$ Radius of adsorbent pellet [m]

$\tau$Dimensionless time

$C_{p,a}$Specific heat capacity of the adsorbed phase [J/(mol.K)]
$C_{p,g}$Specific heat capacity of the gas [J/(mol.K)]
$C_{p,s}$Specific heat capacity of the adsorbent [J/(kg.K)]
$C_{p,w}$Specific heat capacity of the adsorbent [J/(kg.K)]
$k_i$Mass transfer coefficient of component i [s$^{-1}$] or Total number of predictions by SciML
$h_{in}$Inner heat transfer coefficient [W/(m$^2$.K)]
$h_{out}$Outer heat transfer coefficient [W/(m$^2$.K)]
$K_w$ Thermal conductivity of the wall [W/(m.K)]
$K_z$ Effective gas thermal conductivity [W/(m.K)]
$L_{wf}$ Length of column affected by the waterfront [m]
$MW_i$ Molecular weight of component $i$ [kg/mol]
$N$ Number of columns
$r_i$Column inner radius [m]
$r_o $Column outer radius [m]

$t $ Time
$T_a $ Ambient temperature [K]
$T_0 $ Feed temperature [K]
$T_w $ Column wall temperature [K]
$T^*$ Dimensionless temperature
$T^*_{w} $ Dimensionless wall temperature
$U$ Internal energy [J/mol]
$w_l $ Water mole fraction in the feed gas to PSA unit
$y_i $ Mole fraction of component $i$ in the gas phase

$L$Lagrangian function, [-]
$\lambda$Vector of Lagrange multipliers, [-]
$\boldsymbol{\theta}_j^*$Optimal vector of operating conditions, [-]
$\boldsymbol{\theta}_k$Vector of particles, [-]
$n_i$Number of objective functions, [-]
$n_k$Number of particles, [-]
$n_j$Point number on the Pareto front, [-]
$n_{\theta}$Number of decision variables, [-]
$F$Test Fisher, [-]

§ INTRODUCTION

Cyclic adsorption processes are recognized for efficiently separating complex mixtures while presenting energy efficiency and environmental sustainability <cit.>. Among these processes, the pressure swing adsorption (PSA) unit stands out for its diverse range of applications in gas-phase separations. These applications span various industries, including the petrochemical sector <cit.>, $CO_2$ capture <cit.>, and the fine chemicals field <cit.>. Nevertheless, the PSA model exhibits complex dynamics characterized by cyclic behavior instead of steady-state operation. This occurs due to the alternation between the different stages that define the system, which continually vary throughout the cycle. Furthermore, PSA systems often incorporate multiple adsorption beds operating in series or parallel, each with individual characteristics. The quest for energy optimization in pressure swing adsorption (PSA) modeling and the necessity to account for transient effects and multiphase transport presents a challenge. This demands the development of sophisticated approaches that accurately represent the complex behavior of PSA processes and enhance their efficiency. Crucially, these approaches must be designed to optimize process performance while ensuring minimal computational effort, striking a balance between detailed representation and practical feasibility.

Building on the challenges highlighted in the pursuit of energy optimization in PSA modeling, the operation of Pressure Swing Adsorption itself is inherently complex due to the continuous alternation of steps in the process. This complexity is further compounded when trying to achieve efficient operation, as it requires the precise definition and coordination of various operational variables, including pressures, temperatures, gas flow rates, and cycle times. Moreover, PSA processes often involve conflicting performance metrics <cit.>, where optimizing one parameter can adversely affect another <cit.>. This interplay of variables and the trade-offs involved underscore the need for advanced optimization strategies that can navigate these complexities while maintaining process efficiency and low computational burden.  For instance, pursuing higher purity and recovery may lead to trade-offs in the process, making optimizing the PSA process a significant challenge, requiring careful balance among various operational factors to achieve the desired performance and overall efficiency <cit.>.

In this scenario, <cit.> focus on the simulation and optimization of an industrial PSA unit for hydrogen purification. The research involves developing a comprehensive mathematical model that accurately represents the dynamics of a commercial PSA unit. This model is then utilized to optimize key operational parameters, such as cycle time and the steps involved in the pressure swing process. <cit.>  addressed the complexity of optimizing pressure swing adsorption processes. Their work presents a method to optimize a PSA unit for syngas purification using porous aminofunctionalized titanium terephthalate MIL-125. The integrates big data analytics and metaheuristic techniques, aiming to effectively address and manage the uncertainties commonly encountered in the optimization process. This approach highlights the computational challenges of using phenomenological models in PSA process optimization, particularly due to the complexity of the partial differential equation systems.

The difficulties in conducting optimizations based on first-principle models are evident in various fields where rigorous modeling results in complex numerical problems. In such scenarios, the challenge lies in balancing the accuracy of detailed models with the feasibility of computational resources and time constraints.  For PSA processes, <cit.> addressed this issue by proposing an innovative approach to the cost-optimal design of PSA systems. It introduces a model that integrates design elements, such as the number of beds and operational sequences, to minimize the annualized cost of PSA separation systems. The model simplifies the complex dynamics of PSA processes into algebraic equations, balancing capital and operating costs <cit.>. On the other hand, <cit.> introduces a Newton-based approach for determining cyclic steady states in PSA systems, integrating design constraints for more accurate optimization. The study also employs advanced sequential quadratic programming (SQP)-based algorithms for designing optimal PSA processes, incorporating detailed adsorption models. A key innovation is the parallelization of sensitivity calculations, significantly improving computational efficiency.  <cit.> highlight these issues focusing on optimizing a four-step adsorption process for $CO_2$ capture using a novel finite volume-based simulator. This research addresses the numerical challenges in adsorption process simulations, comparing various high-resolution schemes for accuracy and efficiency. The study's core contribution lies in its rigorous optimization approach, which efficiently identifies optimal operating conditions for $CO_2$ capture, balancing process demands with energy consumption. This work marks a significant step in advancing numerical methods for optimizing adsorption processes.

<cit.>, address the computational challenges in simulating and optimizing pressure swing adsorption systems. The authors highlight that optimizing these systems, due to the coupled stiff PDE system and steep fronts moving with time, poses significant computational difficulties for existing optimization techniques. To overcome this, the study develops a reduced-order model (ROM) based on proper orthogonal decomposition (POD). This ROM serves as a low-dimensional approximation to the dynamic PDE-based model, offering a more cost-efficient and manageable approach for the simulation and optimization of PSA systems.

In this scenario, more recently, artificial intelligence (AI), specifically artificial neural networks (ANN), has been proposed as a surrogate model to reduce the computational effort. This approach facilitates the extraction and storage of vast information from databases while requiring minimal computational resources for simulation. Consequently, ANN models have emerged as a bridge connecting rigorous phenomenology with the low computational effort required for real-time applications <cit.>.

For instance,
<cit.> introduced a framework for deep neural networks (DNN) models in PSA units, showcasing the advantages of DNNs in capturing intricate system behavior or representing system performance parameters. Addressing the computational effort related to the optimization procedure of these processes and the issue of properly selecting the optimal sorbent,  <cit.> presents a strategy for simultaneous material screening and process optimization. This strategy leverages the capabilities of deep neural networks to extract knowledge from databases and employs a nested optimization loop that simultaneously considers process and material perspectives. The single objective optimization problem results are analyzed using a Fisher–Snedecor test, which assesses the uncertainties of optimal points and helps build feasible operating regions for the processes.  <cit.> introduces an innovative approach using physics-based artificial neural networks (PANACHE) to simulate and optimize chromatographic processes. This method addresses the computational challenges of solving chromatography's hyperbolic PDEs and nonlinear adsorption isotherms. The approach significantly improves the accuracy of process simulations by learning the underlying PDEs through a physics-constrained loss function. The study demonstrates the effectiveness of this method with binary solute mixtures and various isotherm systems, achieving up to 250 times faster computational speeds.  <cit.> addresses the complexities in simulating and optimizing adsorption-PSA. This study introduces an ANN generalized and adsorbent-agnostic framework, significantly accelerating the simulation and optimization process. The ANN framework facilitates rapid adsorbent screening and process optimization, effectively reducing the computational effort typically associated with these procedures.  <cit.> proposed DNN identification and integration into an economic predictive control scheme, serving as a predictive model within the controller.

In this context, a series of recent studies have delved into applying Deep Neural Network strategies for optimizing  PSA units. Collectively, these works signify the advent of DNN-based approaches in PSA optimization, showcasing the integration of advanced machine learning techniques in refining the efficiency and effectiveness of PSA processes. This trend highlights a  shift towards leveraging data-driven methodologies in the field of PSA optimization <cit.>.

These surrogate models can be useful tools not only for optimization but also for tasks such as online soft-sensing and predictive modeling for model predictive controllers.  In <cit.>, the authors address issues related to PSA complex dynamics. They propose the development of data-driven sensors for the PSA inference problem. The study particularly focuses on addressing the uncertainties inherent in these sensors, proposing a novel approach using soft sensors and deep learning techniques. In <cit.>, the authors extend their DNN-based approach, developing an adaptive digital twin for pressure swing adsorption systems based on neural network models. This research likely explores integrating a novel feedback-tracking system with online learning and uncertainty assessment to improve the performance of PSA systems. In this context, a digital twin suggests a sophisticated simulation model that mirrors the real-world PSA system, allowing for real-time monitoring, prediction, and optimization.

While the aforementioned studies have significantly advanced the literature, a common limitation is the lack of systematic evaluation of optimization optimality. While the advantages of ANNs in these contexts are clear, assessing their accuracy in achieving the optimization tasks is crucial. This is particularly important considering that surrogate-based optimization, like that involving ANNs, can sometimes introduce artificial minima into the optimization landscape. Such occurrences can mislead the optimizer, potentially yielding unreliable results. Therefore, failing to evaluate these models rigorously could lead to representations that do not accurately reflect the system's behavior in practical scenarios, underscoring the need for thorough validation of ANN-based optimization strategies.

Tackling the optimization task while simultaneously identifying a feasible operation region (FOR) remains a complex yet crucial challenge in surrogate-based optimization for PSA processes. This innovative concept, initially introduced by <cit.> for true moving bed reactors, presents a valuable strategy. The FOR approach generates a 'cloud' of optimal solutions, from which a subset can be selected based on specific criteria. These selected points can undergo experimental validation or rigorous model testing, comprehensively evaluating the effectiveness of the optimization. Implementing this methodology offers several benefits, particularly in enhancing understanding the interactions between various operational variables in intricate, multi-stage processes like PSA. Moreover, the FOR concept proves highly beneficial in real-time applications, where the accuracy and reliability of these models are important. This approach not only aids in ensuring the precision of surrogate models but also contributes significantly to the field of PSA process optimization.

In this study, we introduce a novel framework for the multi-objective optimization of pressure swing adsorption units, utilizing deep neural network models and incorporating an assessment of optimality through feasible operation regions. This framework integrates an operations map with the Pareto front, from which samples are selected and rigorously tested against a detailed model to confirm the optimization's effectiveness. This approach significantly reduces the computational effort compared to traditional PDE-based optimization methods while also ensuring the reliability of the surrogate-based optimization results. The key contributions of this work are twofold: firstly, it presents an innovative method for optimizing PSA processes using DNNs, aligned with a Pareto-optimal FOR front; secondly, it provides a comprehensive evaluation of the optimality of DNN-based optimization strategies. Through this study, we aim to enhance the efficiency and accuracy of PSA unit optimization, particularly in the context of carbon capture applications.

§ METHODOLOGY

The first step of the methodology proposed in this study is obtaining the dataset. These datasets are then fed to the surrogate model identification framework. In sequence, a multi-objective optimization proposed is solved using the surrogate models to identify the Pareto Front. Then, the points in the Pareto front are tested to evaluate the optimality of the surrogate-based optimization.

The Fisher-Snedecor test establishes a feasible operational region near the Pareto front. This operational region is subsequently divided into three distinct subgroups, forming the core components for constructing the operational map, as depicted in Figure <ref>. Finally, a detailed process analysis is carried out for the three different regions. Then, the operation map is defined. Further elaboration of these details is provided below.

Diagram illustrating the methodological steps developed in this work.

§.§ Pressure Swing Adsorption to $CO_2$ Capture

There has been a substantial increase in atmospheric $CO_2$ levels, with an emerging scientific consensus indicating that this rise is a primary contributor to recent global warming. This upward trend in global temperatures has been a focal point in the contemporary climate change discussion. In response, continuous effort has been to develop and implement renewable energy technologies to minimize carbon emissions. Despite these advancements, fossil fuels will remain relevant to the energy matrix for some time. Reducing $CO_2$ emissions has become a key focus in this context, with special attention to techniques such as carbon capture and storage (CCS).
For this purpose, a technological alternative that has attracted interest in recent years is separation by cyclic adsorption. Among these techniques, pressure swing adsorption is renowned for its operation that balances cost-efficiency with effectiveness. Particularly, it demonstrates proficiency in facilitating challenging gas phase separations. Operated cyclically, these systems achieve a pseudo-continuous performance where each bed recurrently operates through a series of different steps. The effectiveness of this technology, in terms of purity and recovery capacity, is crucial for optimizing costs and ensuring efficient capture of $CO_2$ emissions.

In this work, PSA system employing the Skarstrom cycle <cit.> was adopted for the separation of $N_2$ and $CO_2$, utilizing zeolite 13X as the adsorbent material. Zeolite 13X is recognized as one of the high-performance adsorbents in the field of gas separation and purification. Characterized by its unique microporous structure and high adsorption capacity, this zeolite stands out in the selective capture of specific molecules, such as $N_2$ and $CO_2$ <cit.>. This process configuration is widely used in PSA due to its simplified operation and long history of industrial use.  This cycle consists of the following stages: (1) pressurization (Pres), (2) adsorption (Ads), (3) heavy reflux (HR), (4) counter-current depressurization (Depres), and (5) light reflux (LR).

In the initial pressurization stage, a gas mixture comprising $CO_2$ and $N_2$ elevates the pressure of the column. In the subsequent adsorption stage, $CO_2$ gets adsorbed in the bed, while $N_2$ is produced. As the bed approaches saturation, an additional high-purity $CO_2$ stream is introduced into the column to purge residual $N_2$ within the bed (HR stage). Following this in the counter-current depressurization stage, the column undergoes depressurization to desorb the  $CO_2$ from the bed, transitioning from a high-pressure state to a low-pressure state. In the light reflux stage, high-purity $N_2$ is injected into the bed to regenerate the column by removing part of the remaining $CO_2$.

In a PSA process, a wide range of operating conditions must be specified throughout the cycles to ensure high purity and recovery of the separated gases. The set of variables may include pressure, feed temperature, step durations, and flow rates, making the optimization of this unit challenging. This set of decision variables must be adjusted to meet operation targets while satisfying purity and recovery constraints.

On the other hand, several performance metrics must be evaluated. The purity ensures that the final product meets the required quality standards.
The recovery metric is important because it determines how much the desired component is recovered from the process, minimizing losses and optimizing efficiency. Due to complex dynamics, obtaining high purity and recovery might not be a simple task.

Illustrative diagram of the pressure swing adsorption steps used in this study

§.§ Phenomenological Model

In this work, the above PSA was modeled by a system of partial differential equations (PDEs) and nonlinear algebraic equations to represent the system's mass, energy, and momentum balances <cit.>. Some assumptions were considered, as shown below:

- Ideal gas law is used to describe the gas phase;

- No concentration, pressure, or temperature gradient in the radial or azimuth directions;

- Linear driving force is used to describe the gas diffusion into the adsorbent;

- Adsorbent properties and void fraction are constant throughout the column;

-  Viscosity of the gas is independent of pressure;

- There is thermal equilibrium between the adsorbent and the gas phase;

- The Ergun equation is used to describe the pressure drop across the bed;

- The column operates adiabatically; wall energy balance is not needed;

- An axially dispersed plug flow model represents bulk fluid flow.

All of the equations employed in the model were converted into dimensionless forms. When working with dimensionless quantities, simulations' numerical stability and accuracy are improved, mitigating problems associated with extremely large or small numbers. The dimensionless variables are provided below:

\begin{align}
&P^*=\frac{P}{P_0},\; T^*=\frac{T}{T_0}, x_{i\ }=\frac{q_i}{q_{s0}}, { \nu^*}_Z=\frac{\nu_Z}{v_0},
\tau=\frac{t\nu_0}{L},
x_i^\ast=\frac{q_i^\ast}{q_{s0}}, Z\ =\frac{z}{L}
\label{adimensionalizacao}
\end{align}

where $P$ is the pressure of the system, $P^*$ is dimensionless pressure, $P_{0}$ is adsorption pressure, $T^*$ is dimensionless temperature, $T$ is temperature, $T_{0}$ is feed temperature, $x_i$ is dimensionless molar loading of component $i$ in the solid phase, $q_i$ is molar loading of component $i$ in the solid phase, $q_{s0}$ is molar loading scaling factor, $v_Z$ is superficial gas velocity, $ \nu^*_Z$ dimensionless superficial gas velocity, $v_0$ is velocity scaling factor, $\tau$ is dimensionless time, $t$ is time, $L$ is the length of the column, $x_i^\ast$ is dimensionless equilibrium molar loading of component $i$ in the solid phase, $q_i^\ast$ is the equilibrium molar loading of component $i$ in the solid phase, $q_{s0}$ is the molar loading scaling factor, $Z$ is the dimensionless length coordinate, $z$ is bed length coordinates.

The mass balance (Equation <ref>) implemented in this model was computed about to the molar fraction of $CO_2$ (represented by $y_i$) in the gas phase. This mass balance is calculated by:

\begin{align}
\frac{\partial y_i}{\partial\tau}=\frac{1}{Pe}\left(\frac{\partial^2y_i}{\partial Z^2}+\frac{1}{P^*}\frac{\partial P^*}{\partial Z}\frac{\partial y_i}{\partial Z}-\frac{1}{T^*}\frac{\partial T^*}{\partial Z}\frac{\partial y_i}{\partial Z}\right)  -{\bar{v}}_Z\frac{\partial y_i}{\partial Z}  +\frac{\Psi T^*}{P^*}\left(\left(y_i-1\right)\frac{\partial x_i}{\partial\tau}+y_i\frac{\partial x_{i+1}}{\partial\tau}\right)
\label{mass_balance}
\end{align}

where $\Psi\ =\left[\frac{(1-\varepsilon)}{\varepsilon}\right]\ \left(\frac{RT_0q_{s0}}{P_0}\right)$, and the Peclet number is calculated by:

\begin{equation}
\label{Peclet}
\end{equation}

The axial dispersion coefficient, $D_L$, is calculated using the following:

\begin{equation}
D_L=0.7D_m+r_p\cdot v_0
\label{eq_DL}
\end{equation}

where $D_m$ is the molecular diffusivity of $CO_2$
- $N_2$ mixture and $r_p$ is the radius of the adsorbent pellet.
The total system pressure was calculated using the bed total mass balance equation:

\begin{align}
\frac{\partial P^*}{\partial\tau}=&\left(-P^*\frac{\partial\left({\bar{v}}_Z\right)}{\partial Z}-{\bar{v}}_Z\frac{\partial\left(P^*\right)}{\partial Z}+{\bar{v}}_Z\frac{P^*}{T^*}\frac{\partial\left(T^*\right)}{\partial Z}\right) \Psi T^*\sum_{i}\frac{\partial\bar{x_i}}{\partial\tau}+\frac{P^*}{T^*}\frac{\partial T^*}{\partial\tau}
\label{balanco_massa}
\end{align}

The mass transfer between the gaseous and solid phases is calculated by the linear driving force (LDF) given by:

\begin{equation}
\frac{\partial x_i}{\partial\tau}=\frac{k_iL}{v_0}(x_i^\ast-x_i)
\label{balanco_solido}
\end{equation}

where $k_i$ is the mass transfer coefficient of component $i$.

The pressure drop throughout the column is calculated using the Ergun equation:

\begin{align}
-\frac{\partial P^*}{\partial Z}=&\frac{150\mu\left(1-\varepsilon\right)^2L}{4r_p^2\varepsilon^2P_0}{\bar{v}}_Zv_0 +\frac{1.75\left(1-\varepsilon\right)L}{2r_p\varepsilon P_0}\left(\sum_{i} y_iMW_i\rho_g\right){{\bar{v}}_Z}^2v_0^2
\label{ergun}
\end{align}

where $\mu$ is gas viscosity, $\varepsilon$ is bed void fraction, $r_p$ is radius of adsorbent pellet and $MW_i$ is molecular weight of component $i$.

The density of the gas phase ($\rho_g$) is given by:

\begin{equation}
\rho_g=\frac{P^*P_H}{RT_0T^*}
\label{massa_esp}
\end{equation}

where $R$ is the universal gas constant and $P_H$ is the high pressure.

The energy balance equation for the gas phase is defined as follows:

\begin{equation}
\frac{\partial T^*}{\partial\tau}=\ \pi_1\frac{\partial^2T^*}{\partial\tau^2}-\ \pi_2{\bar{v}}_Z\frac{\partial T^*}{\partial\tau}+\sum_{i}{(\pi_{3i}+\ \pi_4T^*)\frac{\partial x_i}{\partial\tau}}
\label{balanco energia}
\end{equation}

\begin{equation}
\pi_1=\frac{K_Z}{\left(\rho_gC_{p,g}\varepsilon+(1-\varepsilon)\ \left(C_{p,s}\rho_s+C_{p,a}q_{s0}\right)\right)v_0L}
\label{pi1}
\end{equation}

\begin{equation}
\pi_2=\frac{\rho_gC_{p,g}\varepsilon}{\left(\rho_gC_{p,g}\varepsilon+(1-\varepsilon)\ \left(C_{p,s}\rho_s+C_{p,a}q_{s0}\right)\right)}
\label{p2}
\end{equation}

\begin{equation}
\pi_{3i}=\frac{(1-\varepsilon )(-\Delta H_i)q_{s0}}{T_0 \left ( \rho_g C_{p,g}  \varepsilon + (1-\varepsilon)(C_{p,s} \rho_s + C_{p,a} q_{s,0}) \right )}
\label{p3}
\end{equation}

\begin{equation}
\pi_4=\frac{(1-\varepsilon)(C_{p,g}-C_{p,a})q_{s0}}{\left(\rho_gC_{p,g}\varepsilon+(1-\varepsilon)\ \left(C_{p,s}\rho_s+C_{p,a}q_{s0}\right)\right)}
\label{p4}
\end{equation}

\begin{equation}
\Delta H_i =\Delta U_i-RT^*T_0
\label{Hi}
\end{equation}

where $C_{p,s}$ is specific heat capacity of the adsorbent, $C_{p,g}$ is specific heat capacity of the gas, $C_{p,a}$ is specific heat capacity of the adsorbed phase, $K_Z$ is effective gas thermal conductivity, $\rho_s$ is density of adsorbent, $\Delta H_i$ is enthalpy change, $\Delta U_i$ is internal energy change, $\pi$ is dimensionless groups in the energy balance equations.

This work considers the zeolite 13X to have two temperature-dependent Langmuir sites. We employed isothermal parameters to depict the adsorption process, as summarized in Table <ref>. The associated isotherm equations, derived from <cit.>, are as follows:

\begin{equation}
\label{solido}
\end{equation}

in which

\begin{align}
B_i=b_i exp^{\frac{-\Delta U_{b,i}}{RT}},
\label{parametro_B} \\
D_i=d_i exp^{\frac{-\Delta U_{d,i}}{RT}}.
\label{Parametro_D}
\end{align}

Table <ref> presents some of the parameters and specifications required to solve the phenomenological model for PSA.

Isotherm parameters for $CO_2$ and $N_2$ for zeolite 13X on two sites named b and d.

Parameters                                                                   $CO_2$                                    $N_2$                                    1c     1c
Saturation loading for site b, $q^b_{sat,i}(mol/kg)$       6.21                                   11.9                                  1c     1c

Saturation loading for site d, $q^d_{sat,i}(mol/kg)$       7.15                                   0.00                                  1c     1c

Pre-exponential factor for site b, $b_i (kPa^{-1})$            $4.47\times10^{-8}$     $1.40\times10^{-7}$

Pre-exponential factor for site d, $d_i (kPa^{-1})$          $ 4.71\times10^{-10}$     0.00

Heat of adsorption for site b, $\Delta U_{b,i} (kJ/mol)$     -37.8                                  -19.4

Heat of adsorption for site d, $\Delta U_{d,i} (kJ/mol)$     -37.8                                  0.00

Table <ref> presents the model's initial conditions. As we initiate the simulation, the bed and column are assumed to have reached equilibrium with the gas phase at atmospheric temperature. During the pressurization phase, the same pressure as in the purge phase is maintained. Following the initialization phase, the initial conditions (refer to Equation <ref>) for each subsequent phase are presumed to align with the bed profile after the prior phase. This ensures a transition and coherence of the simulation.

\begin{align}
&P^* = P^*_L,\;    y_i = y_{feed},\;    T^* = T_a,\;   x_i = x_i^*\mid_{y_{feed}}\;
\label{Cond_inicial}
\end{align}

PSA model parameters

Parameters                                                     Value
Length of the column, $L$ $(m)$                                    1

Inlet gas $CO_2$ mole fraction, $y_0$                          0.15

Viscosity of gas, $\mu$ $(Pa.s)$                                 $1.72 \times 10^{-5}$

Void fraction, $\varepsilon$                                   0.37

Molecular diffusivity, $D_m$ $(m^2/s)$          $1.3 \times 10^{-5}$

Thermal conduction in the gas phase,$ k_z$ $(W/m.k)$    0.09

Specific heat of gas, $C_{pg}$ $(J/mol.k)$                       30.7

Specific heat of adsorbed phase,  $C_{pa}$ $(J/mol.k)$           30.7

Radius of the pellets, $r_p$ $(m)$                               $1 \times 10^{-3}$

Specific heat capacity of the adsorbent, $C_{ps}$ $(J/kg.K)$     1070

Molar loading scaling factor, $q_s$ $(mol/kg)$                   5.84

Mass transfer coefficient for $CO_2$, $k_{CO_2}$ $(1/s)$         0.1631

Mass transfer coefficient for $N_2$, $k_{N_2}$ $(1/s)$           0.2044

Purge Pressure, $P_L$ $(Pa)$                                     $1 \times 10^{4}$

Adsorption pressure, $P_0$ $(Pa)$                                $1\times 10^{5}$

Feed temperature of flue gas, $T_0$ $(K)$                       313.15

To solve the system of partial differential equations, we begin by discretizing the spatial variables. In this process, we employ the finite volume method (FVM) and utilize an essentially non-oscillatory scheme known as WENO. Our discretization involves subdividing the equations into 10 volume elements. This value was defined based on a sensitivity analysis. Subsequently, we solve the resulting ordinary differential equations (ODEs) using MATLAB's rigid solver, ode15s. This discretization and equation resolution process ensures an effective and precise solution for the system of PDEs under study, allowing for a detailed analysis of its characteristic dynamic behavior.
This numerical solver was run on a computer featuring an AMD Ryzen 9 7950X processor, 128 GB of RAM, and an RTX 4090 GPU.

PSA model boundary conditions

3lCycle step
1c2*1. Pressurization     1c$Z = 0^+$     1c$P^* = P^*_L$     1c$y_i = y_{feed}$     $T^* = 1$

1c                                       $Z = 1^-$                          $\frac{\partial P^*}{\partial Z} = 0$                         $\frac{\partial{y_i}}{\partial Z} = 0 $              $\frac{\partial T^*}{\partial Z} = 0$

2*2. Adsorption                             1c$Z = 0^+$     1c$P^* = P_high$     1c$y_i = y_{feed}$     $T^* = 1$

$Z = 1^-$                          $P^* = 1$                         $\frac{\partial\bar{y_i}}{\partial Z} = 0 $                         $\frac{\partial T^*}{\partial Z} = 0$

2*3. Heavy reflux                                  1c$Z = 0^+$     1c$P^* = P_high$     1c$y_i = y_{feed}$     $T^* = 1$

$Z = 1^-$                          $P^* = 1$                         $\frac{\partial\bar{y_i}}{\partial Z} = 0 $                         $\frac{\partial T^*}{\partial Z} = 0$

2*4. Depressurization                               1c$Z = 0^+$     $P^* = P^*_L$             $\frac{\partial{y_i}}{\partial Z} = 0 $                        $\frac{\partial T^*}{\partial Z} = 0$

$Z = 1^-$                         $\frac{\partial P^*}{\partial Z} = 0$                          $\frac{\partial{y_i}}{\partial Z} = 0 $                          $\frac{\partial T^*}{\partial Z} = 0$

2*5. Light reflux                                  1c$Z = 0^+$     $P^* = P^*_L$                         $\frac{\partial{y_i}}{\partial Z} = 0 $                          $\frac{\partial T^*}{\partial Z} = 0$

$Z = 1^-$                          $P^* = P^*_L$              $y_i = y_{i,ads}$                          $T^* = T^*_{ads}$

§.§ Data Acquisition and Data Processing

The phenomenological model of the PSA process played an essential role as a virtual plant for generating synthetic data. The controlled generation of synthetic data was done through the design of the experiments (DoE) allowing the exploration of various scenarios and operational conditions, enabling the surrogate model training in challenging situations that may be difficult to encounter in practice. This, in turn, leads to surrogate models that are robust and capable of handling a wider range of real-world scenarios, contributing to improved performance. To do so, we considered the PDE model as a virtual plant accessed through a software-in-the-loop (SIL).

Hence, the first step in obtaining and collecting synthetic data was the design of the experiments. The experimental conditions were defined based on a pseudo-random sampling technique called latin hypercube sampling (LHS) <cit.>. LHS provides a solid statistical foundation for the experiments conducted in the virtual plant. This method combines random and stratified sampling techniques to generate sets of nearly random samples of various variables, avoiding unwanted correlations in the input variables of the DNN model.

The LHS's first step involves sampling each variable's cumulative probability density function. Subsequently, the variable samples are organized based on their correlations to prevent undesired correlations <cit.>.
Therefore, the input variables of the surrogate model (in this case: $t_{press}$, pressurization step time; $t_{depress}$, depressurization step time; $t_{ads}$, adsorption step time; $t_{LH}$, light reflux step time; $t_{HR}$, heavy reflux step time) were defined using the latin hypercube sampling technique to be incorporated into the virtual plant and generate synthetic data regarding $CO_2$ purity ($Pur_{CO_2}$), $CO_2$ recovery ($Rec_{CO_2}$) in the cyclic steady-state (CSS). The minimum limits define the limits for each operational variable
$u_{min}$ = [10, 10, 10, 10, 10]s and maximum limits $u_{max}$ = [1000, 1000, 1000, 1000, 1000]s.

Following the generation of the input matrix using the latin hypercube sampling method, a correlation assessment among the variables was conducted. Figure <ref> illustrates a heat map of all the model's input variables. The results reveal that the correlations exhibit values near zero, suggesting that the input space provided by LHS was thoughtfully designed. This effectively mitigates data distortions and, consequently, ensures that the training process remains free from unwanted biases.
The LHS was used to design 3,000 experiments. Afterward, these disturbances were introduced to the phenomenological model to produce datasets constituting the model's output matrix. Figure <ref> showcases the input representation to introduce disruptions in the phenomenological model. Each experiment was subsequently implemented in the PSA system until it achieved a cyclic steady state.

The data set obtained from the DoE was divided randomly into three subsets according to the following proportion: 75% for training data, 15% for validation data, and 15% for test data.
This work structured the data using a multiple-input single-output (MISO) strategy.

Heatmap of correlations among LHS input signals.

Input representation for LHS-designed experiments for phenomenological model perturbation.

§.§ DNN Identification

After the data collection and preparation, the next step is the configuration of the non-linear function approximator, which in this case is the DNN model. DNN has several parameters to be estimated. One essential step before estimating these parameters is defining the model structure, the hyperparameters. These variables dictate the training and performance of a DNN model. The hyperparameters can be divided into two sets: model and algorithm parameters. The number of layers, neurons, activation function, and layer type are in the first set. In the second set are the tuning parameters of the optimization algorithm, which include the learning rate, momentum, learning rate decay, and several epochs for gradient-based optimizers. All those hyperparameters together form a space of both discrete and continuous variables. Due to their number and different nature, their proper estimation is a complex task. However, optimally adjusting these hyperparameters can enhance the network's generalization capability and mitigate the risk of overfitting.

The usual approach in the machine learning field is trial-and-error-based methods, such as random and grid search for hyperparameter tuning <cit.>. Due to the nature of these methods, they are computationally heavy and not efficient. On the other hand, the most recent in this field was made a few years ago by <cit.>, where a method called HYPERBAND was proposed. This method was demonstrated to be significantly faster in search speed than random grid search while providing a statistical basis for a consistent search of the best set of parameters. This method is based on a formulation of the hyperparameter optimization problem as a pure-exploration problem. It uses a present resource allocated to perform random samples of the search space. The method performs a cuddled loop where an inner loop based on the SUCESSIVEHALVING method <cit.> is coupled to the outer loop based on a grid search. The application of the HYPERBAND requires the definition of the maximum number of epochs and the discarded proportion of non-ideal configurations in each. The HYPERBAND algorithm is one of the most advanced strategies in the field.

Before applying the HYPERBAND algorithm for optimization, the initial step involves defining the hyperparameters of interest and their corresponding search spaces. These hyperparameters may encompass the learning rate, the number of layers in the neural network, batch size, and other factors influencing the model's performance. Each hyperparameter is bounded by a specific search space, typically representing a range of possible values for that parameter. By defining these search spaces, the HYPERBAND algorithm can explore various combinations of hyperparameters to find the optimal set that maximizes the model's performance.
In this study, the initial learning rate, the number of dense layers, the activation functions in each layer, and the number of neurons in each layer were selected as the focus hyperparameters in the quest for the ideal set of parameters for the model. In this study, Table <ref> describes the hyperparameter search space.

Hyperparameter search space for HYPERBAND and best hyperparameters for each performance indicator for DNN.

4cPerformance parameter models

Models    Hyperparameters     Search space     Best hyperparameters

6*[origin=c]90$CO_2$ purity     Initial learning rate     $1\times{10}^{-4}$,$1\times{10}^{-3}$, $1\times{10}^{-2}$     $1\times{10}^{-3}$

Number of dense layers     {1,2,3,4,5,6}     5

Activation function in the dense layers     {relu, tanh}     {relu, relu, relu, relu, linear}

Number of neurons in the dense layer     from 50 to 300 with 10 step     {290, 100, 280, 200, 1}

Number of parameters for the layers     –     {1740, 29100, 28280, 56200, 201}

6*[origin=c]90$CO_2$ recovery     Initial learning rate     $1\times{10}^{-4}$,$1\times{10}^{-3}$, $1\times{10}^{-2}$     $1\times{10}^{-3}$

Number of dense layers     {1,2,3,4,5,6}     4

Activation function in the dense layers     {relu, tanh}     { relu, relu, relu, linear}

Number of neurons in the dense layer     from 50 to 300 with 10 step     {220, 60, 80, 1}

Number of parameters for the layers     –     {1320, 13260, 4880, 81}

The HYPERBAND algorithm iteratively adjusts the model's hyperparameters throughout the optimization process to enhance its performance for the defined objective function. This objective function is evaluated using the training and validation data with the hyperparameters selected within the specified search space. In this context, the function encompasses a loss function that assesses the model's performance on the training and validation datasets.

§.§ DNNs Training and Testing

With the selection of DNN structures, the subsequent phase involves estimating and testing the models' internal parameters. The ADAM algorithm <cit.> was used to train the models based on the configuration given in Table <ref>. Following the training process, a testing step is performed using a set of experiments that were not part of the model training process. The testing results are illustrated in Figure <ref>, enabling a visual assessment of the fitness of each model.
More specifically, the parity plot visually compares model predictions and actual test data. This chart assesses the model's performance and ability to capture system behavior accurately. The random distribution of points along the graph's diagonal line highlights the model's accuracy. The parity plot demonstrates that the AI model's predictions consistently align with the test data across the entire test range.

$CO_2$ Purity

$CO_2$ Recovery

Parity graphics of the model test.

Table <ref> presents the AI model's performance results. To assess the models' performance, we employed standard metrics such as mean absolute error (MAE) and mean squared error (MSE), widely recognized for measuring the discrepancy between actual values and predictions <cit.>. Furthermore, the identification and statistical analysis of a model's predictions can be important in tracking the potential limitations of the model and assessing the need to readjust its parameters.

The results reveal that the models exhibit significantly low values for both MAE and MSE concerning the test data. This attests to the successful model identification and their ability to make accurate predictions. These low error levels underscore the models' reliability and capacity to capture the PSA model's performance metrics. Based on the results in Table <ref>, the obtained error is acceptable and does not justify the increased computational effort.

Final MAE and MSE to test dataset.

1lVariable     MSE                               MAE
$CO_2$ Purity                     $1.13 \cdot 10^{-6}$       $7.44 \cdot 10^{-4}$

$CO_2$ Recovery                   $1.14 \cdot 10^{-6}$     $6.98 \cdot 10^{-4}$

§.§ Multi-objective Optimization Problem for PSA unit

This section presents the mathematical design of the optimization problem proposed here. A multi-objective strategy was followed as $CO_2$ capture in a PSA unit must accommodate the trade-off between unit performance parameters. Instead of exclusively focusing on a single metric, multi-objective optimization considers both objectives simultaneously. This assists in finding solutions that represent the best compromise between purity and recovery, considering the specific constraints and requirements of the process.
Well-adjusted step times can significantly impact the final product's purity and the desired component's recovery. A longer time in the adsorption step may increase purity but potentially reduce recovery. Conversely, a shorter time may result in higher recovery but reduced purity. Therefore, optimizing step times is essential to achieve the desired balance between purity and recovery, tailoring the process to specific application needs. In this work, step times were used as decision variables in the optimization problem.
All process outputs were evaluated in cyclic steady state (CSS) so that the previously mentioned constraints are end-point constraints. The inputs, decision variables $u$, are the same inputs from the data-based model; they were subjected constraints
$[u_{min},u_{max}]$ imposed along the trajectory until the CSS is reached.

In mathematical terms, the optimization problem at hand can be calculated as per the following problem:

\begin{equation}
\begin{cases}
&max \, \, \text{Pur}_{CO_2}(u, \beta),  \\
&max \, \,  \text{Rec}_{CO_2}(u, \beta),
\end{cases}
\label{funcao_obj}
\end{equation}

subject to:

\begin{align}
& u \in [u_{\text{min}}, u_{\text{max}}] \\
& 0\% \leqslant \text{Pur}_{CO_2}(u, \beta), \text{Rec}_{CO_2}(u, \beta) \leqslant 100\%
\label{restricoes}
\end{align}

\begin{equation}
Pur_{CO_2}\ =\frac{{moles}_{CO_2,\left.out\right|Depres}+{moles}_{CO_2,\left.out\right|LR}}{{moles}_{total,\left.out\right|Depres}+{moles}_{total,\left.out\right|LR}}
\label{pureza}
\end{equation}

\begin{equation}
Rec_{CO_2}\ =1-\frac{{moles}_{CO_2,\left.out\right|Ads}}{{moles}_{CO_2,\left.in\right|Pres}+{moles}_{CO_2,\left.in\right|Ads}}
\label{recuperacao}
\end{equation}

Equations <ref> - <ref> are the necessary equations referring to the number of moles at different steps to calculate purity and recovery
of $CO_2$ over the cycles <cit.>.

\begin{equation}
{moles}_{CO_2,\left.out\right|Depres}=\frac{P_0L}{RT_0}\varepsilon A\int_{0}^{\tau_{Depres}}{{\bar{v}}_Zy_{CO_2}\frac{P^*}{T^*}d\tau\ \ \ \ \forall Z=0^+}
\label{moles_1}
\end{equation}

\begin{equation}
{moles}_{CO_2,\left.out\right|LR}=\frac{P_0L}{RT_0}\varepsilon A\int_{0}^{\tau_{LR}}{{\bar{v}}_Zy_{CO_2}\frac{P^*}{T^*}d\tau\ \ \ \ \forall Z=0^+}
\label{moles_2}
\end{equation}

\begin{equation}
{moles}_{total,\left.out\right|Depres}=\frac{P_0L}{RT_0}\varepsilon A\int_{0}^{\tau_{Depres}}{{\bar{v}}_Z\frac{P^*}{T^*}d\tau\ \ \ \ \forall Z=0^+}
\label{moles_3}
\end{equation}

\begin{equation}
{moles}_{total,\left.out\right|LR}=\frac{P_0L}{RT_0}\varepsilon A\int_{0}^{\tau_{LR}}{{\bar{v}}_Z\frac{P^*}{T^*}d\tau\ \ \ \ \forall Z=0^+}
\label{moles_4}
\end{equation}

\begin{equation}
{moles}_{CO_2,\left.out\right|Ads}=\frac{P_0L}{RT_0}\varepsilon A\int_{0}^{\tau_{Ads}}{{\bar{v}}_Zy_{CO_2}\frac{P^*}{T^*}d\tau\ \ \ \ \forall Z=1^-}
\label{moles_5}
\end{equation}

\begin{equation}
{moles}_{CO_2,\left.in\right|Pres}=\frac{P_0L}{RT_0}\varepsilon A\int_{0}^{\tau_{Pres}}{{\bar{v}}_Zy_{CO_2}\frac{P^*}{T^*}d\tau\ \ \ \ \forall Z=0^+}
\label{moles_6}
\end{equation}

\begin{equation}
{moles}_{CO_2,\left.in\right|Ads}=\frac{P_0L}{RT_0}\varepsilon A\int_{0}^{\tau_{Ads}}{{\bar{v}}_Zy_{CO_2}\frac{P^*}{T^*}d\tau\ \ \ \ \forall Z=0^+}
\label{moles_7}
\end{equation}

To solve the optimization problem outlined in Equation <ref> - <ref>, we applied a constrained sliding particle optimization algorithm, as proposed by <cit.>. This algorithm has proven to be a practical choice for optimizing cyclic adsorption processes.

To develop a solid assessment and gain important insights into the decision-making process, the $CO_2$ productivity and the energy required for the CCS cycle were calculated, as shown in the following equations.

\begin{equation}
CO_2 \; productivity = \frac{mol\; CO_2 \; in \; the \; product}{mass \; of \; adsorbent \times cycle \; time}
\label{productivity}
\end{equation}

\begin{equation}
Energy \; required = \frac{energy \; required \; for \; all \; steps }{mass  \; of \; CO_2 \; collected \; in the \; product \; per \; cycle}
\end{equation}

Specifically in terms of the optimization algorithm, this research implemented a variant of PSO, the constraints-based sliding particle swarm optimizer (CSPSO), as proposed by <cit.>, a metaheuristic algorithm grounded on population principles and multilocal minima search. The selection of CSPSO was driven by its ability to handle local minima and to draw the process feasible operating regions. Within the CSPSO algorithm, particles, constituted by decision variables, navigate the problem’s search space. Each particle’s position signifies a potential solution correlated with the set of decision variables. Consequently, PSO utilizes elementary mathematical operations to direct each particle towards areas of promise within the search space. This progression considers the optimal positions found by the individual particle - termed as the individual component, and the most favorable value discovered across all particles, termed as the social component. With each iteration, a new set of positions is obtained. The optimization procedure within PSO persists until a termination condition, such as achieving a maximum number of
iterations or obtaining a satisfactory solution, is met. Upon satisfaction of the termination condition, the algorithm concludes, delivering the optimal solution found during the optimization process.

§.§ Feasible Operation Region Map

In this study, we repurposed the results obtained from assessing the objective function during the optimization problem-solving process to construct a feasible operating region of the process closely aligned with ideal conditions. To build the FOR, we utilized the extended Fisher-Snedecor statistical test. This test was employed in the context of the multi-objective optimization problem addressed in this research to outline the region near the Pareto front.
The Equation <ref> shows the Fisher-Snedecor test applied to this study, with 95% confidence level ($\alpha$) <cit.>.

\begin{equation}
L_i (\boldsymbol{\theta}_k,\boldsymbol{\lambda_k}) \leq L_i (\boldsymbol{\theta}_j^*,\boldsymbol{\lambda_j^*} )+ \frac{(n_k n_j)}{(n_k-n_{\boldsymbol{\theta}}+n_i )}  F_\alpha (n_j,n_k-n_{\boldsymbol{\theta}}+1)
\label{Fisher}
\end{equation}

In this context, the CSPSO optimizer was used with 50 iterations and 500 particles. Figure <ref> illustrates the solution of the optimization problem using the data-driven model. The Pareto front obtained comprises 660 points, presenting significant diversity within the set of solutions.

The optimality of the surrogate-based optimization was evaluated by introducing the decision variables corresponding to the 660 points from the Pareto front in the rigorous model. Hence, if there is an overlap between the rigorous model outputs and the Pareto front found in the surrogate-based optimization, we can consider it a reliable result.

This overlap guarantees the good results achieved in the optimization conducted through the substitute model based on the identified DNN. It was possible to use the optimal points obtained by the DNN model as input to the phenomenological model to establish a comparison with the substitute model. Consulting the phenomenological model takes approximately 600 seconds while using the empirical model requires only 0.1 seconds. This disparity in computational effort between both strategies justifies the option for the substitute approach. On the other hand, we are reducing the phenomenological optimization requirement to an on-demand evaluation of the rigorous model.
For instance, when using 50 particles and conducting 50 iterations, the DNN-based optimization shows convergence in 250 seconds, while the phenomenological model could take 259200 seconds, about 1000X more than the surrogate-based optimization. This notable discrepancy highlights the computational advantage of the DNN model. However, as presented here, this advantage must be used with caution. Therefore, we proposed to sample only the Pareto front of the surrogate-based optimization; this extra step adds 68429 seconds to the optimization procedure, which is still 273 times lower than the phenomenological-based optimization.

Overlap of the Pareto optimal front for the DNN model and phenomenological model and feasible operation region defined by the Fisher-Snedecor test.

With the optimality verified, we can proceed to explore the surrogate-based results. In this work, we propose dividing the FOR into three groups using clustering techniques. Each subgroup will represent a zone on the Pareto front where each conflicting objective is favored. region 1 prioritizes one objective, region 3 prioritizes another, whereas region 2 represents a balance between the two. Furthermore, the operation region clusters provide flexibility for the process operation, thus constituting a tool capable of assisting in decision-making based on the Pareto map clusters. This tool can play an important role in process monitoring, providing relevant information about the behavior of the system's optimal conditions.

Therefore, the clustering method aggregates the population points within the FOR according to a metric. Among various possible metrics, the weighted average distance (WPGMA) was used, in which the (Euclidean) distance between the clusters is calculated as a simple average. The method aims to construct subgroups of data based on defined centroids, which will form each of the three regions. This is done by calculating the distance between the centroid and the surrounding points. The points that compose each region are those with the smallest distance from the corresponding centroid.

The next section will delve into the insights, comparisons, and knowledge gained regarding the system through the process operation map.

§ DISCUSSION

One of the objectives of this work is to use the operation regions defined based on the methodology presented in this study to discuss the process behavior in the context of the resulting operating map.
Using the operational region derived from the Pareto front facilitates the process operating in terms of guiding the resource allocation. Hence, these results are produced as a decision-making support tool.

In the optimization problem presented, there are two objectives, and the objective function is continuous. From the decision-maker's perspective, there are three options: an action that prioritizes one objective over the other, the opposite action, or an action that would strike a balance between the two objectives. Thus, to represent these possible courses of action, three clusters were used.

The FOR was subdivided into three distinct operational regions presented in Figure <ref>, each with a specific focus and priorities. The first of these zones, denoted as region 1, primarily emphasizes maximizing $CO_2$ purity, striving to achieve the highest possible levels of purity within the constraints of the process. In region 2, the objective shifts to striking a delicate balance between $CO_2$ purity and $CO_2$ recovery, where both factors are optimized concurrently to meet specific operational requirements. Finally, region 3 prioritizes maximizing $CO_2$ recovery without considering the $CO_2$ purity. The division of FOR into these operational zones offers a more nuanced approach to process optimization, allowing for tailored solutions that align with the specific objectives and constraints of the system. Hence, the results are converted into a guideline for decision-making.

Feasible operation regions defined by clustering.

The operation map is also conducted from the perspective of decision variables associated with each respective cluster. This aspect is crucial as it allows for a more comprehensive evaluation of how the process can be effectively operationalized to achieve the desired objectives.

Figure <ref> depicts the combinations of decision variables, their interrelationships, and how regions 1, 2, and 3 are associated with the operating variables.

It is important to note that in region 1, where $CO_2$ purity is the primary focus, the duration of the light reflux step is limited to around 100 seconds. This implies the process’s sensitivity to longer LR times, which could introduce impurities into the light reflux output stream and, consequently, reduce the purity of $CO_2$. In regions 2 and 3, the duration of the light reflux step increases, leading to a gradual decrease in $CO_2$ purity.

Furthermore, the interplay between LR time and $CO_2$ purity significantly affects process efficiency. In region 1, the LR time duration given by our methodology ensures high $CO_2$ purity but may limit the overall $CO_2$ recovery. In contrast, Regions 2 and 3 allow longer LR times, prioritizing higher $CO_2$ recovery but at the expense of slightly lower $CO_2$ purity. This delicate balance between purity and recovery is a crucial factor in optimizing the process and tailoring it to specific application requirements.

Assessing $CO_2$ purity from the standpoint of the adsorption step, it is evident that higher $CO_2$ purities require extended durations in this phase to ensure efficient $CO_2$ capture and reduce the presence of $N_2$ in the bed. Thus, during the depressurization and light reflux steps, the stream is enriched in $CO_2$.
On the other hand, $CO_2$ recovery is affected by the extended duration of the adsorption step. The longer $CO_2$ contaminates the output stream of the adsorption step, the lower the $CO_2$ recovery will be. These relationships become evident when revisiting Equations <ref> and <ref>.
The balance between purity and recovery is a critical factor in optimizing the PSA process. Achieving the desired purity while maximizing $CO_2$ recovery requires carefully adjusting the cycle durations.

For the remaining Press, Depres, and HR steps, a wide range of combinations is possible within a broad time frame, ranging from 10 to 1000 seconds. This suggests a sloppiness in these variables regarding the optimal conditions.

On the other hand, the variability within this wide time range offers process engineers the flexibility to adapt the PSA process to different operating conditions and specific objectives. However, the emphasis on the correlation and constraints of $t_{LR}$ and $t_{ads}$ underscores the pivotal role of these parameters in determining the overall performance and trade-offs in the PSA process.

Feasible operation map from the point of view of decision variables.

A three-dimensional (3D) analysis, as represented in Figure <ref>, offers another perspective when evaluating decision variables. This methodology allows observation of particle dispersion, highlighting the interconnection of three operational variables within distinct clusters. The duration of the process steps appears to be a key factor that influences the purity and recovery of $CO_2$.

In this context, a sharp contrast is observed in the dispersion of data points in different regions. Region 1, which prioritizes $CO_2$ purity, presents less dispersion for the three variables. In contrast, region 3, which emphasizes the recovery of $CO_2$, shows a significantly higher degree of dispersion, suggesting a wider range of variable interactions.
This distinction highlights the specificity and precision required in adjusting operational variables to achieve high purity levels.

Feasible operation map in 3D for $t_{depres}$, $t_{ads}$ and $t_{LR}$, (region 1 - yellow, region 2 - blue, region 3 - red)

To enrich our analysis, we present in Figure <ref> a representation of the system's dynamics in the CSS. As shown in Table <ref>, we selected a point in each operation region within the operations map, allowing us to evaluate the dynamic behavior at the end of each step. Table <ref> shows the corresponding performances for each point in the operating regions.

This operational profile suggests that while region 1 prioritizes $CO_2$ purity over recovery with moderate energy consumption, region 2 prioritizes $CO_2$ recovery with slightly better productivity. region 3 maximizes $CO_2$ recovery at the expense of higher energy demand and lower purity. The selection of the optimal operation region will depend on the specific objectives of the process and cost-benefit considerations related to energy consumption and the desired product quality.

Operational parameters     Region 1     Region 2     Region 3
Time of pressurization step / s     503.6     663.1     1001

Time of depressurization step / s     367.5     465.7     488.8

Time of adsorption step / s     989.6     1001.0     635.6

Time of light reflux step / s     59.7     160.1     215.2

Time of heavy reflux step / s     526.8     521.0     558.8
$CO_2$ purity / %     80.0     73.1     51.8

$CO_2$ recovery / %     73.6     89.9     99.4

Energy consumption / kWh     125.5     210.7     320.4

$CO_2$ productivity / (mol/kg.h)     0.338     0.358     0.251

Operation parameters and their corresponding performance for each representative point for the regions of the operation map.

The duration of each stage for each selected point is highlighted in Figure <ref>. This figure represents the concentration of the column output at the end of a cycle in the CSS.
In Figure <ref>, it is noteworthy that region 1 exhibits significantly higher $CO_2$ concentrations than regions 2 and 3. Specifically, in the depressurization and light reflux stages, the $CO_2$ molar fraction is notably higher when compared to the points in regions 2 and 3, leading to higher $CO_2$ purity. In Figure <ref>, we can also observe that the adsorption step exhibits a higher fraction of $CO_2$ in the column outlet. This indicates that the $N_2$-enriched stream is being contaminated, compromising $CO_2$ purity. Note that at the Pareto front, purity does not reach values higher than 80%.

History of $CO_2$ mole fractions for a cycle in the CSS for each selected point in each operating region. (region 1 - yellow, region 2 - blue, region 3 - red).

Figure <ref> displays the molar fraction profiles of $CO_2$ and $N_2$ along the bed at the end of each stage in the CSS. The aim is to compare these profiles at three distinct points representing each region in the operational map. In region 1, contamination in the enriched $N_2$ flow during adsorption becomes evident. Notably, the $CO_2$ fraction does not reach zero at the upper end of the bed, indicating the persistence of residual $CO_2$ in the enriched $N_2$ stream at the stage exit. This contamination diminishes the efficiency of $CO_2$ recovery. Another important factor is related to region 3, where the recovery of $CO_2$ is prioritized. It is observed that the capture of $CO_2$ occurs only in 50% of the bed, while the remaining part of the bed adsorbs the lighter component, indicating an underutilization of the column. In region 2, 90% of the bed is dedicated to capturing $CO_2$, while 10% captures $N_2$.

From the standpoint of $CO_2$ purity, at the point corresponding to region 1, there was a more excellent capture of $CO_2$ compared to the other points. This results in a purer $CO_2$ stream at the exit of depressurization, favoring the overall cycle's purity.

Internal bed profiles for  mole fraction of $CO_2$ and $N_2$ (solid line $CO_2$ fraction, dashed line $N_2$ fraction).

The depressurization and light reflux stages are crucial for desorbing the heavy products and regenerating the column adsorption capacity. This is achieved by reducing the pressure in the system, thus, taking advantage of the adsorption isotherm.  The energy consumed in the desorption-related stage primarily depends on the necessity for a vacuum and for how long it is rewuired.

Thus, increasing recovery in a PSA system can lead to a higher energy demand during the depressurization and purge stages. This extension in the adsorption stage typically implies higher operating times in a high-pressure condition to ensure more efficient retention.

Figure <ref> shows that the points representing regions 2 and 3 depict longer depressurization and light reflux stages than those corresponding to region 1. This pattern also holds for the pressurization and adsorption stages. When the durations of these stages are combined, they yield a higher total time for regions 2 and 3. As a result, there has been an increase in energy requirements, as illustrated in Figure <ref>.

In Figure <ref>, we highlight the energies required for the depressurization and light reflux stages, which stand out as the most energy-intensive aspects. The energies required for the other stages are comparatively less significant.

Energy required per cycle in the CSS and for the depressurization and light reflux stages for the points representing regions 1, 2, and 3.

Another essential discussion concerns the inherent trade-off between purity, recovery, and productivity. To achieve a high purity of $CO_2$, the PSA system is operated to retain $CO_2$ and eliminate $N_2$ selectively. However, this pursuit of purity often results in lower productivity, as a significant portion of the gas feed is discarded to ensure high purity.
On the other hand, when productivity is the focus, the PSA system is optimized to capture as much $CO_2$ as possible without necessarily achieving very high purity. The other side of the trade-off is the fact that recovery refers to the percentage of $CO_2$ captured relative to the total $CO_2$ available in the source. At the same time, productivity concerns the total amount of $CO_2$ captured per cycle time in the CSS (Equation <ref>).
In general, achieving the desired recovery requires the system to be more efficient in $CO_2$ capture. This often involves the application of stricter operating conditions that ensure a higher purity of $CO_2$ in the final product. However, an increase in recovery often comes at the cost of reduced productivity due to the extended time required for each complete cycle in PSA.
Conversely, when productivity is the primary priority, the system is optimized to capture as much $CO_2$ as possible, even if this means lower recovery.

In Figure <ref>, productivity is depicted for each point corresponding to regions 1, 2, and 3. Observing this graph makes it striking how it illustrates the trade-off between purity, recovery, and productivity. Regions that require higher purity and recovery often display lower productivity. In contrast, region 2 stands out for achieving a balance between purity and recovery, resulting in a more even productivity. This graph underscores the importance of finding the ideal combination of these factors, considering the specific process goals, to achieve the desired performance in the $CO_2$ capture system.

Comparison between the energy required and $CO_2$ productivity in the CSS for each representative point in the operating region.

In Figure  <ref>, a comparison is also made between the energy required for the three regions and productivity. It is noted that increasing productivity often involves the pursuit of more efficient operating conditions, such as shorter cycle times. However, these improvements often demand more energy, as longer gas compression, adsorbent regeneration, and other process steps may require additional energy. Therefore, an increase in productivity can lead to an increase in energy consumption.

On the other hand, reducing energy consumption is an important goal in many industries, as it lowers operational costs and environmental impacts. However, energy reduction is often associated with decreased productivity since operating processes under milder or less efficient conditions may produce fewer desired components within a given time frame.

Thus, higher $CO_2$ purity (region 1) will require less energy consumption and yield higher productivity, whereas greater recovery (region 3) entails higher energy requirements and lower productivity.

§ CONCLUSIONS

This study introduces a methodology of surrogated optimization of cyclic adsorption units while offering a robust strategy for assessing the effectiveness of the optimization process. This methodology was evaluated in practice to optimize a pressure swing adsorption unit for $CO_2$ capture.

To do so, we identified four DNN models within a multiple-input, single-output framework. Each of these models was capable of predicting critical process performance parameters. Subsequently, these models were integrated into an optimization framework, utilizing a particle swarm optimizer and a statistical evaluation of the resulting population. The outcome was a comprehensive representation of the Pareto front, including its regions of feasible operation.

To validate the effectiveness of the optimization process, decision variables corresponding to the Pareto front were tested in the phenomenological model, confirming the reliability of the surrogate model. The overlap ensures successful optimization through the substitute model based on the identified DNN, highlighting a significant computational advantage. While the phenomenological model requires substantial time for certain scenarios, the DNN achieves convergence much more efficiently, underscoring the effectiveness of the surrogate-based approach.

Finally, we analyzed the feasible operational domains of the decision variables, accompanied by a correlation map detailing their interplay. This conclusive analysis shed light on the most influential variables shaping process behavior and provided insights into how they can be manipulated to attain specific optimal criteria.

The outcomes of this study offer a significant contribution to decision-making processes. It presents a comprehensive operational map that guides operators in determining the optimal process location and prioritizing specific objectives. This strategic tool enhances the ability to make informed choices, effectively aligning operations with targeted goals.

§ DECLARATION OF COMPETING INTEREST

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

§ ACKNOWLEDGEMENT

The present work contributes to the completion of a sub-project at SUBPRO, a research-based innovation center within Subsea Production and Processing at the Norwegian University of Science and Technology. The authors would like to express their gratitude for the financial support received from SUBPRO, funded by the Research Council of Norway through grant number 237893, major industry partners, and NTNU.


# A complete Bernstein function related to the fractal dimension of Pascal’s
pyramid modulo a prime

Christian Berg

(12.2.2024)

###### Abstract

Let $f_{r}(x)=\log(1+rx)/\log(1+x)$ for $x>0$. We prove that $f_{r}$ is a
complete Bernstein function for $0\leq r\leq 1$ and a Stieltjes function for
$1\leq r$. This answers a conjecture of David Bradley that $f_{r}$ is a
Bernstein function when $0\leq r\leq 1$.

Dedicated to the memory of Bent Fuglede (1925-2023)

Mathematics Subject Classification: Primary 30E20; Secondary 26A48

Keywords. Bernstein function, complete Bernstein function, Stieltjes function,
Nevanlinna-Pick function.

##  1 Introduction and main results

For $0<r$ we consider the function

$f_{r}(x):=\frac{\log(1+rx)}{\log(1+x)},\quad x>0.$ (1)

It extends by continuity to the interval $(\max(-1,-1/r),0]$ with the value
$f_{r}(0)=r$. Furthermore, $\lim_{x\to\infty}f_{r}(x)=1$.

In a study of the fractal dimension of certain sets in number theory coming
from the number of multinomial coefficients congruent to a given residue
modulo a prime, David Bradley is led to the conjecture, that $f_{r}$ is a
Bernstein function when $0<r<1$, see [4]. We prove that this is true as part
of the following result:

######  Theorem 1.

The function $f_{r}$ is a complete Bernstein function when $0<r<1$ and a
Stieltjes function when $r>1$.

Bernstein functions are by definition $C^{\infty}$-functions
$f:(0,\infty)\to[0,\infty)$ satisfying

$(-1)^{n-1}f^{(n)}(x)\geq 0,\quad x>0,n=1,2,\ldots.$ (2)

In particular, Bernstein functions are non-negative, increasing and concave.
To verify the inequalities (2) for $f_{r}$ seems a complicated task, because
one has to differentiate a quotient. We prove instead the stronger statement,
that $f_{r}$ belongs to the subclass $\mathcal{CBF}$ of complete Bernstein
functions studied in [8]. Our proof uses complex analysis.

Bernstein functions have a long history, and they play an important role in
probability theory and potential theory, see [8] and [1].

We also need the class ${\mathcal{S}}_{\lambda}$ of generalized Stieltjes
functions of order $\lambda>0$ defined as the functions
$f:(0,\infty)\to[0,\infty)$ having a representation of the form

$f(x)=c+\int_{0}^{\infty}\frac{d\mu(t)}{(x+t)^{\lambda}},$ (3)

where $c\geq 0$, and $\mu$ is a positive measure on $[0,\infty)$ such that the
integral in (3) converges for $x>0$. When $\lambda=1$ the functions in
$\mathcal{S}=\mathcal{S}_{1}$ are just called Stieltjes functions. A function
in the class $\mathcal{S}_{\lambda}$ clearly extends to a holomorphic function
in the cut plane $\mathbb{C}\setminus(-\infty,0]$. For more information on
these functions see [7]. We mention in particular that
$\mathcal{S}\subset\mathcal{S}_{2}$.

Complete Bernstein functions are characterized by six equivalent conditions
given in Theorem 6.2 of [8]. We shall establish condition (iv), stating that
they are Pick functions with certain properties. For the convenience of the
reader we have included some key results about Pick functions in Section 2.

We denote by $\operatorname{Log}$ the principal branch of the logarithm, which
is holomorphic in the cut plane $\mathbb{C}\setminus(-\infty,0]$ and defined
by

$\operatorname{Log}(z)=\log|z|+i\operatorname{Arg}(z),\quad\operatorname{Arg}(z)\in(-\pi,\pi).$

The expression

$f_{r}(z):=\frac{\operatorname{Log}(1+rz)}{\operatorname{Log}(1+z)}$

is a holomorphic extension of $f_{r}$ to the cut plane
$\mathbb{C}\setminus(-\infty,\max(-1,-1/r)]$, ($z=0$ being a removable
singularity).

We first give more details about the case $0<r<1$.

######  Theorem 2.

For each $0<r<1$ the function $f_{r}$ is a complete Bernstein function given
for $z\in\mathbb{C}\setminus(-\infty,-1]$ by

$f_{r}(z)=\frac{\operatorname{Log}(1+rz)}{\operatorname{Log}(1+z)}=r+\int_{1}^{\infty}\frac{z}{z+t}\sigma_{r}(t)\,dt,$
(4)

where

$\sigma_{r}(t)=\left\\{\begin{array}[]{cl}\frac{-\log(1-rt)}{t\left(\log^{2}(t-1)+\pi^{2}\right)},&1<t<1/r,\\\
\infty,&t=1/r,\\\
\frac{\log((t-1)/(rt-1))}{t\left(\log^{2}(t-1)+\pi^{2}\right)},&1/r<t.\end{array}\right.$

On the interval $(-1,\infty)$ the function $f_{r}$ is strictly concave and
strictly increasing from $0$ to $1$. The function $1-f_{r}(z)$ is a Stieltjes
function given by

$1-f_{r}(z)=\int_{1}^{\infty}\frac{t\sigma_{r}(t)}{z+t}\,dt,$ (5)

and $f_{r}(z)/z$ is a Stieltjes function given by

$f_{r}(z)/z=\int_{0}^{\infty}\frac{d\mu_{r}(t)}{z+t},\quad\mu_{r}=r\delta_{0}+1_{(1,\infty)}(t)\sigma_{r}(t)\,dt,$
(6)

where $\delta_{0}$ is the Dirac measure at $0$. The function $\sigma_{r}$
tends to 0 for $t\to 1^{+}$ and for $t\to\infty$. The following formulas hold

$\int_{1}^{\infty}\sigma_{r}(t)\,dt=1-r,\quad\int_{1}^{\infty}\frac{\sigma_{r}(t)\,dt}{1+t^{2}}=\frac{\log
2\log(1+r^{2})+\pi\arctan(r)}{\log^{2}(2)+\pi^{2}/4}-r.$ (7)

We include a Maple plot of $\sigma_{r}$.

The function $\sigma_{r}(t),1<t<1/r,\,1/r<t<5$, where $r=1/2$.

######  Remark 3.

There exists a number $r_{0}\approx 0.1$ such that $\sigma_{r}(t)$ is
increasing on $(1,1/r)$ when $r_{0}\leq r<1$, but when $0<r<r_{0}$ then
$\sigma_{r}(t)$ is decreasing in a certain subinterval
$[a_{r},b_{r}]\subset(1,1/r)$ and otherwise increasing.

We next give more details about the case $r>1$.

######  Theorem 4.

For $r>1$ the function $f_{r}$ is a Stieltjes function given for
$z\in\mathbb{C}\setminus(-\infty,-1/r]$ by

$f_{r}(z)=1+\int_{1/r}^{\infty}\frac{\omega_{r}(t)}{z+t}\,dt,$ (8)

where

$\omega_{r}(t)=\left\\{\begin{array}[]{cl}\frac{-1}{\log(1-t)},&1/r\leq
t<1,\\\ 0,&t=1,\\\
\frac{\log((rt-1)/(t-1))}{\log^{2}(t-1)+\pi^{2}},&1<t.\end{array}\right.$

On the interval $(-1/r,\infty)$ the function $f_{r}$ is strictly convex and
strictly decreasing from $\infty$ to $1$ and

$\int_{1/r}^{\infty}\frac{\omega_{r}(t)}{t}\,dt=r-1.$ (9)

We include a Maple plot of $\omega_{r}$.

The function $\omega_{r}(t),1/r\leq t<1,\,1<t<1.5$, where $r=2.5$.

######  Remark 5.

If $f(z)$ is a complete Bernstein function or a Stieltjes function, then so is
$f(cz)$ for any $c>0$. In particular, for $0<a<b$

$f_{a/b}(bz)=\frac{\operatorname{Log}(1+az)}{\operatorname{Log}(1+bz)}$

is a complete Bernstein function, and

$f_{b/a}(az)=\frac{\operatorname{Log}(1+bz)}{\operatorname{Log}(1+az)}$

is a Stieltjes function.

Proofs of the theorems are given in Section 3.

##  2 Preliminaries about Pick functions

Let $\mathbb{H}:=\\{z\in\mathbb{C}\mid\Im z>0\\}$ denote the open upper half-
plane.

A holomorphic function $f:\mathbb{H}\to\mathbb{C}$ is called a Pick function
if $\Im f(z)\geq 0$ for all $z\in\mathbb{H}$.

For a Pick function $f$ the imaginary part $\Im f$ is a non-negative harmonic
function, and it has the mean value property. It is therefore either
identically zero or strictly positive on $\mathbb{H}$. This means that a Pick
function is either a real constant function or it maps $\mathbb{H}$ into
itself. Pick functions are treated in [5], [8] and in the Appendix to [9].

A Pick function $f$ has the representation

$f(z)=\alpha+\beta
z+\int_{-\infty}^{\infty}\left(\frac{1}{t-z}-\frac{t}{1+t^{2}}\right)\,d\mu(t),\quad
z\in\mathbb{H},$ (10)

where $\alpha\in\mathbb{R},\beta\geq 0$ and $\mu$ is a positive measure on
$\mathbb{R}$ satisfying the integrability condition
$\int(1+t^{2})^{-1}\,d\mu(t)<\infty$. The triple $(\alpha,\beta,\mu)$ is
uniquely determined by $f$ as

$\alpha=\Re f(i),\quad\beta=\lim_{y\to\infty}f(iy)/iy,\quad\lim_{y\to
0^{+}}(1/\pi)\Im f(x+iy)\,dx=\mu,$ (11)

where the last convergence is in the vague topology for measures. In fact,
using convolution we can write

$(1/\pi)\Im f(x+iy)=(1/\pi)\beta y+P_{y}\star\mu(x),\quad
P_{y}(t)=\frac{1}{\pi}\frac{y}{t^{2}+y^{2}},$

and it is well-known that the Poisson kernel $P_{y}$ converges to the Dirac
measure $\delta_{0}$ when $y\to 0$.

This implies in particular, that if there is an open interval
$I\subseteq\mathbb{R}$ such that

$\lim_{y\to 0^{+}}(1/\pi)\Im f(x+iy)=\varphi(x),\quad x\in I$ (12)

and the convergence is uniform on compact subsets of $I$, then the restriction
of $\mu$ to the interval $I$ has the continuous density $\varphi$ with respect
to Lebesgue measure.

The measure $\mu$ has mass $m\geq 0$ at a point $a\in\mathbb{R}$ if and only
if

$\lim_{y\to 0^{+}}y\Im f(a+iy)=m.$ (13)

In fact,

$y\Im f(a+iy)=\beta
y^{2}+\int\frac{y^{2}}{(t-a)^{2}+y^{2}}\,d\mu(t)\to\mu(\\{a\\})$

for $y\to 0^{+}$ because

$\frac{y^{2}}{(t-a)^{2}+y^{2}}\leq\frac{2}{1+(t-a)^{2}},\quad
0<y<1,t\in\mathbb{R},$

and therefore $2(1+(t-a)^{2})^{-1}$ is an integrable majorant with respect to
$\mu$.

For a Stieltjes function $f$ given by (3) with $\lambda=1$, we get

$-(1/\pi)\Im f(-x+iy)=P_{y}\star\mu(x),$

and we get similar results for Stieltjes functions as for Pick functions.

##  3 Proofs

Proof of Theorem 2. We need the following boundary minimum principle for
harmonic functions, see e.g. [6, p.5].

######  Proposition 6.

Let $u$ be a harmonic function in $\mathbb{H}$ satisfying

(i) $\liminf_{z\to t,z\in\mathbb{H}}u(z)\geq 0$ for each $t\in\mathbb{R}$,

(ii) $\liminf_{|z|\to\infty,z\in\mathbb{H}}u(z)\geq 0$.

Then $u(z)>0$ for all $z\in\mathbb{H}$ unless $u$ is identically zero.

From the expression

$f_{r}(z)=\frac{\log|1+rz|+i\operatorname{Arg}(1+rz)}{\log|1+z|+i\operatorname{Arg}(1+z)},\quad
z\in\mathbb{H},$

we see that $f_{r}(z)\to 1$ and hence $\Im f_{r}(z)\to 0$ for $|z|\to\infty$
in $\mathbb{H}$.

The harmonic function $\Im f_{r}(z)$ in $\mathbb{H}$ is given by

$\Im
f_{r}(z)=\frac{\log|1+z|\operatorname{Arg}(1+rz)-\log|1+rz|\operatorname{Arg}(1+z)}{\log^{2}|1+z|+\operatorname{Arg}^{2}(1+z)}.$
(14)

It does not seem possible to see directly that this expression is positive,
but it will follow from Proposition 6. We claim that for $t\in\mathbb{R}$

$\varphi_{r}(t):=\lim_{z\to t,z\in\mathbb{H}}(1/\pi)\Im
f_{r}(z)=\left\\{\begin{array}[]{cl}0,&t\geq-1,\\\
\frac{-\log(1+rt)}{\log^{2}(-t-1)+\pi^{2}},&-1/r<t<-1,\\\ \infty,&t=-1/r,\\\
\frac{\log((-t-1)/(-1-rt))}{\log^{2}(-t-1)+\pi^{2}},&t<-1/r.\end{array}\right.$

For this we use that $\operatorname{Arg}(z)\in(0,\pi)$ for $z\in\mathbb{H}$
and

$\lim_{z\to
t,z\in\mathbb{H}}\operatorname{Arg}(z)=\left\\{\begin{array}[]{cl}0,&t>0,\\\
\pi,&t<0,\end{array}\right.\quad\lim_{z\to 0,z\in\mathbb{H}}\log|z|=-\infty.$

The convergence to $\varphi_{r}(t)$ is uniform for $t$ in compact subsets of
the open intervals $(-\infty,-1/r),(-1/r,-1),(-1,\infty)$, and the expressions
of $\varphi_{r}(t)$ to the right are non-negative for all $t\in\mathbb{R}$.

We have shown that the conditions of Proposition 6 are satisfied for $\Im
f_{r}$, and the $\liminf$ condition is a true limit in all cases. Therefore
$\Im f_{r}(z)>0$ for all $z\in\mathbb{H}$. This shows that $f_{r}$ is a Pick
function with a representation

$f_{r}(z)=\Re
f_{r}(i)+\int_{-\infty}^{\infty}\left(\frac{1}{t-z}-\frac{t}{1+t^{2}}\right)\,d\mu_{r}(t),$
(15)

see Section 2, because

$\beta_{r}:=\lim_{y\to\infty}\frac{f_{r}(iy)}{iy}=0,$

as is easily seen. We have

$\alpha_{r}:=\Re f_{r}(i)=\frac{\log
2\log(1+r^{2})+\pi\arctan(r)}{\log^{2}(2)+\pi^{2}/4}>0.$

The measure $\mu_{r}$ has the density $\varphi_{r}$ in the open intervals
$(-\infty,-1/r),(-1/r,-1)$ and $(-1,\infty)$ because of (12), and there are no
masses at the points $-1/r,-1$ by (13). Replacing $t$ by $-t$ we finally get

$f_{r}(z)=\alpha_{r}+\int_{1}^{\infty}\left(\frac{t}{1+t^{2}}-\frac{1}{t+z}\right)\varphi_{r}(-t)\,dt.$
(16)

In the notation of formula (6.8) in [8] we have
$\varphi_{r}(-t)\,dt=(1+t^{2})d\rho(t)$, and (4) follows with
$\sigma_{r}(t)=\varphi_{r}(-t)/t$. Putting $z=0$ in (16) we get

$r=\alpha_{r}-\int_{1}^{\infty}\frac{\varphi_{r}(-t)}{t(1+t^{2})}\,dt=\alpha_{r}-\int_{1}^{\infty}\frac{\sigma_{r}(t)}{1+t^{2}}\,dt,$

and the second formula in (7) follows.

The formula in Theorem 2 for $\sigma_{r}$ shows that it is integrable over the
interval $(1,\infty)$. From formula (4) we then get

$f_{r}(z)=r+\int_{1}^{\infty}\sigma_{r}(t)\,dt-\int_{1}^{\infty}\frac{t\sigma_{r}(t)}{z+t}\,dt,$

so for $z\to\infty$ on the real axis we get

$1=r+\int_{1}^{\infty}\sigma_{r}(t)\,dt,$

and (5) follows together with the first formula in (7).

From (5) it is clear that $f_{r}$ is strictly increasing and concave on
$(-1,\infty)$. Formula (6) follows immediately from (4). $\square$

######  Corollary 7.

For $0<r<1$ we have

$0<f_{r}(x)<r\frac{1+x}{1+rx},\quad x>-1,$ (17)

and

$\log|1+z|\operatorname{Arg}(1+rz)>\log|1+rz|\operatorname{Arg}(1+z),\quad
z\in\mathbb{H}.$ (18)

###### Proof.

The derivative of $f_{r}$ satisfies

$f_{r}^{\prime}(x)\log^{2}(1+x)=\log(1+x)\frac{r}{1+rx}-\log(1+rx)\frac{1}{1+x},$
(19)

and since $f_{r}$ is strictly increasing, we get (17).

Formula (18) follows from (14) because $\Im f_{r}(z)>0$ for $z\in\mathbb{H}$.
∎

######  Remark 8.

There is an important subclass of $\mathcal{CBF}$ consisting of the Thorin-
Bernstein functions, denoted $\mathcal{TBF}$ and studied in Chapter 8 of [8].
According to Theorem 8.2 of [8] we can conclude that
$f_{r}\notin\mathcal{TBF}$ for $0<r<1$ because a necessary and sufficient
condition for $f_{r}$ to be a Thorin-Bernstein functions is that
$t\sigma_{r}(t)$ is an increasing function for $t>0$, and this is never true,
since the function is infinity for $t=1/r$ and finite for larger $t$ and tends
to 0 for $t\to\infty$.

A tedious calculation shows that $t\sigma_{r}(t)$ is increasing on $(1,1/r)$.
Compare this with Remark 3.

We also recall that $f\in\mathcal{CBF}$ belongs to $\mathcal{TBF}$ if and only
if $f^{\prime}\in\mathcal{S}$.

From (4) or (5) we get

$f_{r}^{\prime}(z)=\int_{1}^{\infty}\frac{t\sigma_{r}(t)}{(z+t)^{2}}\,dt,$

showing that $f_{r}^{\prime}$ is a generalized Stieltjes function of order 2,
which is not surprising as this holds for all complete Bernstein functions. In
conclusion we have $f_{r}^{\prime}\in\mathcal{S}_{2}\setminus\mathcal{S}$ for
$0<r<1$.

We recall that generalized Stieltjes functions of order 2 are logarithmically
completely monotonic by a deep theorem of Steutel and Kristiansen, see Theorem
2.1 in [2] for details. Consequently, $f_{r}$ belongs to the class of Horn-
Bernstein functions, introduced in [3].

The Bernstein representation for $f_{r}$, cf. Equation (6.1) in [8], is given
as

$f_{r}(x)=r+\int_{0}^{\infty}(1-e^{-xs})m_{r}(s)\,ds,\quad x>0,$

where $m_{r}$ is the completely monotonic function

$m_{r}(s)=\int_{1}^{\infty}e^{-ts}t\sigma_{r}(t)\,dt,\quad s>0$

according to Theorem 6.2 in [8].

Proof of Theorem 4. We first recall from Chapter 7 in [8] that for
$f\in\mathcal{CBF}$ which is not identically zero, then $z/f(z)$ also belongs
to $\mathcal{CBF}$ and is called the conjugate of $f$. Furthermore, by Theorem
6.2 in [8] we then get that $1/f\in\mathcal{S}$.

Assume now $r>1$. From (22) we have

$f_{r}(x)=\frac{1}{f_{1/r}(rx)},\quad x>0,$

so $f_{r}$ is the reciprocal of a complete Bernstein function and therefore
$f_{r}\in\mathcal{S}$. Since $f_{r}(x)\to 1$ for $x\to\infty$ we have

$f_{r}(z)=1+\int_{0}^{\infty}\frac{d\nu_{r}(t)}{z+t},\quad z\in\mathbb{H}$

for a positive measure $\nu_{r}$ on $[0,\infty)$. From Section 2 we know that
the positive measures

$-\frac{1}{\pi}\Im f_{r}(-t+iy)\,dt$

converge vaguely to $\nu_{r}$ for $y\to 0^{+}$. From (14) we get for $y>0$

$\displaystyle-\frac{1}{\pi}\Im f_{r}(-t+iy)=$
$\displaystyle\frac{1}{\pi}\frac{\log|1-rt+iry|\operatorname{Arg}(1-t+iy)-\log|1-t+iy|\operatorname{Arg}(1-rt+iry)}{\log^{2}|1-t+iy|+\operatorname{Arg}^{2}(1-t+iy)},$

which for $y\to 0^{+}$ converges to

$\omega_{r}(t)=\left\\{\begin{array}[]{cl}0,&t<1/r,\\\
\frac{-1}{\log(1-t)},&1/r<t<1,\\\ 0,&t=1,\\\
\frac{\log((rt-1)/(t-1))}{\log^{2}(t-1)+\pi^{2}},&1<t,\end{array}\right.$

and the convergence is uniform for $t$ in compact subsets of the open
intervals $(-\infty,1/r),(1/r,1),(1,\infty)$. Using the results of Section 2,
it follows that $\nu_{r}$ has the density $\omega_{r}(t)$ with respect to
Lebesgue measure, and there are no masses at $t=1/r,1$, hence (8) follows.

From (8) we see that $f_{r}$ is strictly decreasing and convex on
$(-1/r,\infty)$, and for $z=0$ we get (9). $\square$

######  Corollary 9.

For $r>1$ we have

$f_{r}(x)>r\frac{1+x}{1+rx},\quad x>-1/r,$ (20)

and

$\log|1+z|\operatorname{Arg}(1+rz)<\log|1+rz|\operatorname{Arg}(1+z),\quad
z\in\mathbb{H}.$ (21)

###### Proof.

Equation (20) follows from (19) using that $f_{r}$ is strictly decreasing.

Equation (21) follows from (14) because $\Im f_{r}(z)<0$ for $z\in\mathbb{H}$,
since $f_{r}$ is a Stieltjes function. ∎

##  4 A convolution equation

The Laplace transform of a positive measure $\mu$ on $[0,\infty)$ is the
function

$\mathcal{L}\mu(x)=\int_{0}^{\infty}e^{-xt}\,d\mu(t),\quad x>0,$

provided that the integral is finite for all $x>0$. If $\mu$ has a continuous
density $f:(0,\infty)\to[0,\infty)$ with respect to Lebesgue measure $m$ on
$(0,\infty)$, i.e., $d\mu=fdm$, we write $\mathcal{L}f=\mathcal{L}\mu$.

With this in mind we have

$\mathcal{L}(\mathcal{L}\mu)(y)=\int_{0}^{\infty}e^{-yx}\mathcal{L}\mu(x)\,dx=\int_{0}^{\infty}\frac{d\mu(t)}{y+t},\quad
y>0.$

Note that the family of functions $f_{r},r>0$ satisfies

$f_{r}(x)f_{1/r}(rx)=1,\quad r,x>0.$ (22)

Assume now that $0<r<1$ and rewrite (22) as

$\frac{f_{r}(x)}{x}f_{1/r}(rx)=\frac{1}{x}=\mathcal{L}m(x)=\mathcal{L}1(x),\quad
x>0.$ (23)

Using the Laplace transformation, Equation (6) can be written

$\frac{f_{r}(x)}{x}=\mathcal{L}(\mathcal{L}\mu_{r})(x)=\mathcal{L}g_{r}(x),\quad
g_{r}(s)=r+\mathcal{L}\sigma_{r}(s),$ (24)

and Equation (8) can be written

$f_{1/r}(rx)=1+\int_{1}^{\infty}\frac{\omega_{1/r}(rt)}{x+t}\,dt=1+\mathcal{L}h_{r}(x),$
(25)

where

$h_{r}(s)=\mathcal{L}\omega_{1/r}(rt)(s)=\int_{1}^{\infty}e^{-st}\omega_{1/r}(rt)\,dt,\quad
s>0.$

Using well-known properties of Laplace transformation, we see that (23) is
equivalent to the convolution equation

$g_{r}\star(\delta_{0}+h_{r})(x)=1,\quad x>0$

or

$g_{r}(x)+\int_{0}^{x}g_{r}(s)h_{r}(x-s)\,ds=1,\quad x>0.$

Acknowledgement The author wants to thank David Bradley for having suggested
the title of the present work, which is inspired by [4]. The author also wants
to thank René Schilling for useful comments to the manuscript.

## References

  * [1] C. Berg and G. Forst, Potential Theory on Locally Compact Abelian Groups, Ergebnisse der Mathematik und ihrer Grenzgebiete Band 87, Springer-Verlag, Berlin-Heidelberg-New York, 1975\.
  * [2] C. Berg, S. Koumandos and H. L. Pedersen, Nielsen’s beta function and some infinitely divisible distributions, _Math. Nachrichten_ , 294 (2021), 426–449. (DOI:10.1002/mana.201900217).
  * [3] C. Berg and H. L. Pedersen, A Family of Horn-Bernstein Functions, _Experimental Mathematics_ , 32:3 (2023), 505–513. (DOI:10-1080/10586458.2021.1980460).
  * [4] D. M. Bradley, On the number of multinomial coefficients congruent to a given residue modulo a prime. Preprint.
  * [5] W. F. Donoghue, Jr., Monotone Matrix Functions and Analytic Continuation. Springer Verlag, Berlin, Heidelberg, New York 1974.
  * [6] J. L. Doob, Classical Potential Theory and Its Probabilistic Counterpart. Grundlehren der mathematischen Wissenschaften 262, Springer-Verlag, Berlin, Heidelberg, New York, 1984.
  * [7] S. Koumandos and H .L. Pedersen, On Asymptotic Expansions of Generalized Stieltjes Functions, _Comput. Methods Funct. Theory_ 15 (2015), 93–115. (DOI: 10.1007/s40315-014-0094-7).
  * [8] R. L. Schilling, R. Song and Z. Vondraček, Bernstein functions. Theory and applications. De Gruyter Studies in Mathematics 37, Second Edition, de Gruyter, Berlin 2012.
  * [9] K. Schmüdgen, The Moment Problem, Graduate Texts in Mathematics Vol. 277. Springer International Publishing AG 2017.

Christian Berg
Department of Mathematical Sciences, University of Copenhagen
Universitetsparken 5, DK-2100 Copenhagen, Denmark
email<EMAIL_ADDRESS>