# Optimal Design of Membrane Cascades for Gaseous and Liquid Mixtures via
MINLP

Jose Adrian Chavez Velasco ‡ Davidson School of Chemical Engineering, Purdue
University, West Lafayette, IN-47906  Radhakrishna Tumbalam Gooty ‡ Davidson
School of Chemical Engineering, Purdue University, West Lafayette, IN-47906
Mohit Tawarmalani† Krannert School of Management, Purdue University, West
Lafayette, IN-47906  Rakesh Agrawal∗ Davidson School of Chemical Engineering,
Purdue University, West Lafayette, IN-47906

###### Abstract

Given the growing concern of reducing CO2 emissions, it is desirable to
identify, for a given separation carried out through a membrane cascade, the
optimum design that yields the lowest energy consumption. Nevertheless,
designing a membrane cascade is challenging since, there are often multiple
feasible configurations that differ in their energy consumption and cost. In
this work, we develop a Mixed Integer Non-linear Program (MINLP) that, for a
given binary separation, which may be either liquid or gaseous, finds the
cascade and its operating conditions that minimize energy consumption. To
model the separation at each membrane in the cascade, we utilize the
analytical solution of a system of differential and algebraic equations
derived from the crossflow model and the solution-diffusion theory. We provide
numerical evidence which shows that our single-stage membrane model accurately
predicts experimental data. Unfortunately, the resulting membrane model is
non-convex and, even state-of-the-art solvers struggle to prove global
optimality of the cascades and the operating conditions identified. In this
paper, we derive various cuts that help with relaxation quality and,
consequently, accelerate convergence of branch-and-bound based solvers. More
specifically, we demonstrate, on various examples, that our cuts help branch-
and-bound solvers converge within 5% optimality gap in a reasonable amount of
time and such a tolerance level was not achieved by a simple formulation of
the membrane model. The proposed optimization model is an easy-to-use tool for
practitioners and researchers to design energy efficient membrane cascades.

Keywords: Energy efficiency, Separations, Membrane cascades

Corresponding authors<EMAIL_ADDRESS><EMAIL_ADDRESS>
‡ These authors contributed equally

##  1 Introduction

Membrane technology has gained considerable interest in recent years.
Distinctive features such as modular scale-up flexibility, operational
simplicity and relatively low capital cost make them attractive for the
separation of a variety of mixtures. A few applications that have seen
commercial success include water desalination at large scale, lube oil
dewaxing, reduction of CO2 in natural gas, distributed scale nitrogen
production from air, and so forth [1, 2]. Besides these examples, a spectrum
of applications exists where membranes can potentially be used on a large
scale (see [1, 2]). The development of advanced materials with enhanced
properties, such as higher area-to-volume ratio, has contributed substantially
towards expanding the potential use of membranes. Despite the progress,
membranes still have at most moderate permeabilities and selectivity values in
a majority of applications. Consequently, multiple membrane modules connected
in a sequence, referred to as a membrane cascade or simply a cascade
hereafter, is needed for separations which require one of the components to be
recovered at high purity and high recovery. Nevertheless, identifying the most
attractive cascade, even for binary separations, remains a challenging task.

It is well-known that, for a given separation, multiple membrane cascades
exist. These cascades differ in terms of the overall energy input and cost [3,
4, 5]. This energy requirement and cost also depend on the operating
conditions of the cascade. The resulting energy consumption of separation is
often a large fraction of the overall energy input of the chemical plant and
this energy is usually sourced from the combustion of fossil fuels.
Consequently, the separation may contribute significantly to the carbon
footprint of the plant. Given rising concerns about greenhouse emissions, it
is essential that energy efficient cascades be designed. Not surprisingly,
this problem has been studied extensively in the literature. Nevertheless, we
show that the current literature has not addressed global optimization of
high-fidelity models. This is particularly important because most models for
the permeation process use nonlinear nonconvex equations, and the local
optimization techniques and/or meta-heuristics used in the literature do not
guarantee that the discovered cascades are globally optimal, or even close in
energy consumption to such a cascade. Here, we provide evidence that sometimes
the local optimization solvers, or simply local solvers hereafter, converge to
a suboptimal solution that requires much higher energy than the global optimal
solution. Besides, our computational experience suggests that, without a good
starting solution, the local solvers fail to converge in a majority of cases.
In addition to finding the most efficient membrane cascade, the global
optimization solution is also important for comparing the membrane technology
against an alternate separation technology. Although global optimization
techniques are available for some separation technologies, such as
distillation [6, 7, 8, 9], they are not yet available for design of membrane
cascades making it difficult to perform head-to-head comparisons. Besides,
since membrane material research is an active research endeavor, global
optimization techniques are needed to ascertain selectivity and permeabilities
of components at which membrane technologies may out-compete other separation
processes, thereby helping to identify separations for which new membrane
cascades are likely to yield most promising results. These reasons motivate
the development of an effective global optimization approach for designing
membrane cascades.

Towards the goal, we present several modeling advances in this article. The
highlights of the article are summarized in the following. First, we describe
a unified permeator model that is applicable for both gaseous and liquid
mixtures. In its default form, the permeator model is a differential-algebraic
equation (DAE) system. We express it as a system of algebraic equations by
solving the DAE system analytically. This is in contrast to the common
practice in the literature, where the DAE system is solved approximately using
various discretization techniques. The analytical solution addresses the
trade-off between the accuracy and the complexity of the model, which we shall
discuss in detail in section 3 (§3). Second, we present a novel unified mixed-
integer nonlinear program (MINLP) that is formulated to identify the membrane
cascade which minimizes the overall power input. We solve the MINLP using an
off-the-shelf global optimization solver, or simply a global solver hereafter,
such as BARON [10, 11]. A global solver has built-in heuristics that generate
multiple good-quality initial points for local search. Further, it continues
local search with new initial points as it explores various parts of the
search region until a proof of optimality is obtained. Once a global solver
converges, the incumbent solution is guaranteed to be within a specified
tolerance from the global optimum. Third, for all of our test cases, global
solvers fail to solve the MINLP to the desired optimality tolerance in a
reasonable amount of time. To address this challenge, we derive additional
cuts for the problem using physical insights and by exploiting the
mathematical properties of the governing equations. Through numerical
examples, we demonstrate that the additional cuts expedite the convergence
characteristics of BARON and we are able to solve the MINLP within
5%-optimality tolerance. Fourth, we apply our methods to two industrially
important applications: separation of propylene/propane mixture and the
separation of p-xylene from a mixture of xylene isomers.

The rest of the paper is organized as follows. We present a survey of the
existing literature in  §2. We describe the unified permeator model for a
single stage in  §3, and present our MINLP formulation along with
computational experiments in  §4. We discuss two case studies in  §5, and
present concluding remarks in  §6.

##  2 Literature review

There are two aspects to designing a cascade for a given separation. First,
identify the candidate set; which we define as the set of all potentially
attractive cascades for a given separation. Second, from the candidate set,
identify the cascade along with its operating condition that optimizes the
desired objective. The approaches used to generate the candidate set can be
broadly classified as follows: (i) heuristic/intuition/empirical observation
based approach, and (ii) superstructure based approach. In the first approach,
the candidate set is identified based on physical insights and engineering
judgements (see [12, 13, 14, 15, 16, 17, 5, 18, 19]). In the second approach,
a superstructure is postulated, and the candidate set is obtained by
discarding appropriate units and/or connections (nodes and/or arcs) from the
superstructure (see [20, 4, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32]).
We will adopt the second approach in this work.

Next, a permeator model is needed to assess a membrane cascade. Marriott and
Sørensen[33] used detailed mass, momentum, and energy balances to model the
transport process on both retentate and permeate sides. Their model also takes
into account concentration polarization and non-ideal solution properties. The
resulting model is a system of partial differential algebraic equations. While
this model accurately describes the process, it is computationally challenging
to use it within an optimization framework. As such, most models in the
literature make various assumptions that simplify the equations describing the
transport process on both sides of the membrane. For example, the crossflow
model of Weller and Steiner [34, 35], when combined with flux equations
derived from the solution-diffusion theory, yields such a model that is used
in the literature. Here, the derivation assumes that the bulk phase
concentration, velocity, and temperature gradients in the transverse direction
(the direction perpendicular to the surface of the membrane) are negligible on
both sides of the membrane. Consequently, the permeator model reduces to a
simple ordinary differential algebraic equation (DAE) system. We use this
simplified permeator model in our current work. Despite the simplifying
assumptions, it has been shown that the model agrees well with experimental
membrane performance data [36, 37]. In this work, we provide additional
evidence that validates the cross-flow model against experimental data. One
reason for a good agreement of the cross-flow model is that most commercial
membranes are either asymmetric membranes or a membrane with a thin dense
layer on a porous support. The works, [38, 39, 40] further accounted for
pressure drop along the flow direction. However, as we will show in §3, even
without this change, the model agrees well with experimental data. Therefore,
we will neglect the effect of pressure drop along the flow direction. Aliaga
et al. [22] and Adi et al. [30] additionally assume that the retentate and
permeate sides are perfectly mixed. This leads to the well-known perfect
mixing model [34], which can be described using a system of algebraic
equations. However, as we will show in §3, the predictions from the perfect
mixing model are less reliable than those from the crossflow model. Besides
assuming perfect mixing, Aliaga et al. [22] also assume that the separation
factor (defined as
$(y^{\textnormal{per}}/(1-y^{\textnormal{per}}))/(x^{\textnormal{in}}/(1-x^{\textnormal{in}})$,
where $y^{\textnormal{per}}$ and $x^{\textnormal{in}}$ denote the mole
fraction of the most permeable component in the permeate and feed streams,
respectively) is constant. One consequence of this assumption is that the
permeate composition depends only on the feed composition. This assumption is
therefore too strong, because the permeate composition depends significantly
on the trans-membrane pressure ratio and the stage cut (fraction of the feed
permeated through the membrane).

Once the candidate set and a permeator model have been decided, one of the
following approaches is used to identify the optimal cascade and its operating
conditions (1) explicit enumeration, (2) meta-heuristic approaches such as
genetic algorithms or simulated annealing, and (3) mathematical programming
approaches.

In the first approach, for each cascade, the degrees of freedom are
identified, and an exhaustive sensitivity analysis is performed over the
admissible operating range of degrees of freedom. The operating condition
which optimizes the desired objective is taken to be its optimal operating
condition. Next, the cascade which maximizes/minimizes the desired objective
at its optimal operating condition is taken to be the optimal cascade (see
[14, 13]). The explicit enumeration makes this approach computationally
expensive particularly when there are many degrees of freedom for each cascade
and the number of cascades is large; for example, cascades containing four or
more stages.

In the second approach, simulated annealing or a genetic algorithm is used to
determine the optimal cascade and its optimal operating condition (see [25,
26, 32]). This approach does not provide a global optimality certificate.

In the third approach, a mathematical program is formulated and solved using
standard solvers (e.g. BARON [10, 11], IPOPT [41], DICOPT [42], etc.). This
approach has been used in [23, 24, 22, 25, 26, 30, 31, 27], each of which
formulates and solves an MINLP to simultaneously determine the optimal cascade
and its optimal operating condition. These formulations use binary variables
to model the presence/absence of arcs connecting different membrane stages. In
contrast, Pathare and Agrawal [21] enumerated all cascades in the candidate
set explicitly and then formulated a nonlinear program (NLP) for each cascade,
which was solved to determine its optimal operating condition. When each
cascade in the candidate set had been optimized, the cascade which minimizes
the desired objective was taken to be the optimal cascade. The mathematical
programming approach suffers from the following challenge. Sensitivity
analysis and meta-heuristic approaches can use black-box models, so the first
two methods can use higher fidelity models of transport phenomena. On the
other hand, mathematical programming approaches typically require explicit
functions and/or gradients, and most solvers available today do not allow DAE
systems as model constraints. This led prior approaches to resort to the use
of discretization procedures to convert the DAE system to a system of
algebraic equations. For example, Uppaluri et al. [25, 26] and Scholz et al.
[27] approximate each membrane stage with several perfect mixers connected in
series. Qi and Henson [23, 24] used a combination of Gauss quadrature and a
fourth order Runge-Kutta-Gill method to convert the DAE system into a system
of algebraic equations. Kunde and Kienle [31] discretized each membrane into
elements, and in each element, they assumed that the molar flux of each
component is linear. This enabled them to compute the permeate and retentate
compositions algebraically. All discretization schemes exhibit a trade-off
between the model complexity (number of nonlinear equations, etc.) and the
accuracy of the solution. The finer the discretization scheme the more
accurate is the representation of the DAE system, but the number of nonlinear
nonconvex equations is also larger. This makes it challenging to obtain an
optimality certificate. Most of the works in the literature use a local solver
to solve their formulation. For instance, Qi and Henson [23, 24] use DICOPT++
[42], Aliaga et al. [22] use SBB [43], Pathare and Agrawal [21] use fmincon in
MATLAB [44]. As we mentioned in §1, these local solvers can get trapped in
suboptimal solutions. We address the above challenges by using the analytical
solution of the DAE system, and by using a global solver, such as BARON. We
express the analytical solution in a different form than it is known in the
literature [34, 35], which facilitates the global solution of the optimization
problem with off-the-shelf global solvers. We remark that BARON was used in
some studies such as those of Adi et al. [30], Scholz et al. [27], and Kunde
and Kienle [31]. However, we note that there are limitations. First, as
mentioned before, the perfect mixing model used by Adi et al. [30] does not
agree well with experimental data, and Scholz et al. [27] do not report
whether optimality certificate was obtained. The exception is Kunde and Kienle
[31] who do report that they obtain an optimality certificate. Nevertheless,
we remark that they treat the trans-membrane pressure difference as a
parameter and do not report extensive numerical results with their model.
Therefore, the scaling and robustness of their solution procedure is not well
understood.

Finally, Agrawal and Xu [15, 16, 45], and Pathare and Agrawal [21]
investigated the effect of the exergy loss due to mixing, referred to as
mixing losses, on the overall compression power. Here, the authors sought the
operating condition that minimizes (or completely eliminates in the case of
[21]) the total mixing loss. Empirically, they observed that this approach
yields a near optimal solution for cascades containing five or more stages.
However, for cascades containing fewer stages, their approach does not always
yield the global optimum.

##  3 Permeator model

In the rest of the article, we denote the most and the least permeable
components as $A$ and $B$, respectively. Further, we denote the separation of
a binary mixture $AB$ as $A$/$B$. We use the solution-diffusion theory [46] to
model the local flux of each component through the membrane, and the crossflow
model, proposed by Weller and Steiner [34, 35], to model the overall
permeation process.

For both gaseous and liquid mixtures, the solution-diffusion theory is the
widely-accepted mechanism of mass transfer through dense polymeric membranes
[46]. According to this theory, the constituent components of the mixture are
separated due to their differences in solubility and diffusivity within the
membrane. Depending on the type of the mixture, the local flux of each
component through the membrane can be obtained using the equations below [46].

For liquids, $\displaystyle\left\\{\begin{aligned}
&n_{A}=\mathit{PM}_{A}\left[x-y\exp\left(-\frac{V_{A}(P^{\textnormal{out}}-P^{\textnormal{per}})}{RT}\right)\right]\\\
&n_{B}=\mathit{PM}_{B}\left[(1-x)-(1-y)\exp\left(-\frac{V_{B}(P^{\textnormal{out}}-P^{\textnormal{per}})}{RT}\right)\right]\end{aligned}\right.,$
(1a) For gases, $\displaystyle\left\\{\begin{aligned}
&n_{A}=\mathit{PM}_{A}\left[P^{\textnormal{out}}x-P^{\textnormal{per}}y\right]\\\
&n_{B}=\mathit{PM}_{B}\left[P^{\textnormal{out}}(1-x)-P^{\textnormal{per}}(1-y)\right]\end{aligned}\right.,$
(1b)

where (i) $n_{A}$, $\mathit{PM}_{A}$ and $V_{A}$ (resp. $n_{B}$,
$\mathit{PM}_{B}$, and $V_{B}$) correspond to the local flux, permeance and
molar volume of $A$ (resp. $B$), (ii) $x$ and $P^{\textnormal{out}}$ (resp.
$y$ and $P^{\textnormal{per}}$) denote the local mole fraction of $A$ and the
total pressure on the retentate (resp. permeate) side, and (iii) $R$ and $T$
denote the universal gas constant and the absolute temperature of the mixture
respectively. The expression inside the square brackets in both (1a) and (1b)
corresponds to the driving force, and it is different for liquids and gases.
We derive a unified expression for the local flux of each component that is
applicable for both gases and liquids. This enables us to formulate a common
optimization model for both gases and liquids. Towards this, we define the
following variables:

$\displaystyle u=\begin{cases}\ln{r},&\text{for a gaseous mixture}\\\ \Delta
P^{trans},&\text{for a liquid mixture}\end{cases},$ (2a)
$\displaystyle\beta=\begin{cases}1,&\text{for a gaseous mixture}\\\
0,&\text{for a liquid mixture}\end{cases},$ (2b) $\displaystyle
C_{A}=\begin{cases}1,&\text{for a gaseous mixture}\\\ V_{A}/RT,&\text{for a
liquid mixture}\end{cases},$ (2c) $\displaystyle
C_{B}=\begin{cases}1,&\text{for a gaseous mixture}\\\ V_{B}/RT,&\text{for a
liquid mixture}\end{cases},$ (2d)

where $r=P^{\textnormal{out}}/P^{\textnormal{per}}$ is the pressure ratio, and
$\Delta P^{trans}=P^{\textnormal{out}}-P^{\textnormal{per}}$ is the trans-
membrane pressure difference. We now obtain a expression for the local fluxes
by expressing (1a) and (1b) in terms of $u$, $\beta$, $C_{A}$ and $C_{B}$ as

$\displaystyle
n_{A}=\mathit{PM}_{A}(P^{\textnormal{out}})^{\beta}\left[x-y\,\mathrm{e}^{(-C_{A}\,u)}\right],$
(3a) $\displaystyle
n_{B}=\mathit{PM}_{B}(P^{\textnormal{out}})^{\beta}\left[(1-x)-(1-y)\,\mathrm{e}^{(-C_{B}\,u)}\right].$
(3b)

Figure 1 shows a schematic of the permeation process across a membrane
employing the cross-flow pattern. We make the following assumptions:

  1. 1.

The pressure drop along the membrane module due to the bulk flow of both
permeate and retentate streams is not substantial.

  2. 2.

Concentration polarization does not occur near the surface of the membrane.

  3. 3.

Mass transfer resistance in the bulk permeate and retentate streams is
negligible.

  4. 4.

The separation takes place isothermally.

  5. 5.

Membrane selectivity is independent of the operating pressure and the
composition of the mixture.

Note that, when the Assumption 1 is valid, the degree of separation obtained
from an asymmetric membrane employing either the counter-current (see Figure
2(a)) or the co-current (see Figure 2(b)) flow pattern is close to that
obtained with the cross-flow pattern [36]. This is because, the porous layer
prevents axial mixing of the local permeate just outside of the dense layer.
Consequently, the flux profiles and the net separation remains the same
regardless of the flow pattern. Therefore, our models are also applicable for
cascades employing either the co-current or the counter-current flow pattern,
provided that the Assumption 1 holds and they use asymmetric membranes.

Figure 1: Conceptual representation of crossflow pattern Figure 2: Permeation
across an asymmetrical membrane under counter-current and co-current flow
patterns. The local composition outside of the dense layer ($y$) does not
depend on the flow pattern.

Let the flowrate, pressure and composition of the mixture entering the
membrane module be $f^{\textnormal{in}}$, $P^{\textnormal{out}}$, and
$x^{\textnormal{in}}$, respectively. Consider a differential control volume as
shown in Figure 1 (dashed rectangle). Mass balance of component $A$ across the
control volume is given by $d(fx)=y\,df$, where $x$ (resp. $y$) corresponds to
the local mole fraction of component $A$ on the retentate (resp. permeate)
side, and $f$ corresponds to the local molar flowrate on the retentate side.
The mass balance equation can be rearranged as

$\frac{dx}{df}=\frac{y-x}{f},\quad
x(f^{\textnormal{in}})=x^{\textnormal{in}},\quad
f\in[f^{\textnormal{in}},f^{\textnormal{out}}],$ (4)

where $x(f^{\textnormal{in}})=x^{\textnormal{in}}$ is the initial condition,
and $f^{\textnormal{out}}$ denotes the net flowrate of the retentate leaving
the membrane module. Since there is no axial mixing in the porous layer, the
local mole fraction of each component on the permeate side is simply the ratio
of the local flux of the component to the total local flux i.e.,
$y=n_{A}/(n_{A}+n_{B})$ and $(1-y)=n_{B}/(n_{A}+n_{B})$. Since both the
equations are linearly dependent, we use only the former. We substitute
$n_{A}$ and $n_{B}$ from (3) and rearrange (see §A.1 for details) to obtain

$\displaystyle y-x=k+k(S-1)y-\frac{k\,S}{S-(S-1)y},$ (5)

where $S=\mathit{PM}_{A}/\mathit{PM}_{B}$ is the selectivity of component $A$
w.r.t component $B$ and

$\displaystyle k=\frac{(S-1)-(Se^{-C_{A}u}-e^{-C_{B}u})}{(S-1)^{2}}.$ (6)

We solve the DAE system in (4) and (5) analytically (see §A.2 for derivation)
to obtain

$\displaystyle
S\ln\frac{y^{\textnormal{out}}}{y^{\textnormal{in}}}-\ln\frac{1-y^{\textnormal{out}}}{1-y^{\textnormal{in}}}-k(S-1)^{2}\ln\frac{y^{\textnormal{out}}-x^{\textnormal{out}}}{y^{\textnormal{in}}-x^{\textnormal{in}}}=k(S-1)^{2}\ln(1-\theta),$
(7)

where $y^{\textnormal{in}}=y(f^{\textnormal{in}})$,
$y^{\textnormal{out}}=y(f^{\textnormal{out}})$,
$x^{\textnormal{out}}=x(f^{\textnormal{out}})$, and stage cut
$\theta=(f^{\textnormal{in}}-f^{\textnormal{out}})/f^{\textnormal{in}}$. Mole
fractions $y^{\textnormal{in}}$ and $y^{\textnormal{out}}$ are related to
$x^{\textnormal{in}}$ and $x^{\textnormal{out}}$ via (5) i.e.,

$\displaystyle
y^{\textnormal{in}}-x^{\textnormal{in}}=k+k(S-1)y^{\textnormal{in}}-\frac{k\,S}{S-(S-1)y^{\textnormal{in}}},$
(8) $\displaystyle
y^{\textnormal{out}}-x^{\textnormal{out}}=k+k(S-1)y^{\textnormal{out}}-\frac{k\,S}{S-(S-1)y^{\textnormal{out}}}.$
(9)

Therefore, given $f^{\textnormal{in}}$, $x^{\textnormal{in}}$, $r$ (or $\Delta
P^{trans}$), and stage cut $\theta$, the mole fraction of component $A$ in the
retentate, $x^{\textnormal{out}}$, can be determined by solving (7)–(9)
simultaneously. Next, the mole fraction of component $A$ in the permeate,
$y^{\textnormal{per}}$, can be obtained from the overall component mass
balance around the membrane module,

$f^{\textnormal{in}}x^{\textnormal{in}}=f^{\textnormal{out}}x^{\textnormal{out}}+f^{\textnormal{per}}y^{\textnormal{per}},$
(10)

where
$f^{\textnormal{per}}=f^{\textnormal{in}}\theta=f^{\textnormal{in}}-f^{\textnormal{out}}$.

###  3.1 Validation of the permeator model

We now show the validity of the model by comparing the predicted permeate and
retentate mole fractions as a function of stage cut against the experimental
data for O2/N2 and CO2/CH4 separations in Figures 3 and 4, respectively. In
the Figures, we also show the predictions obtained from the perfect mixing
model [34]. Clearly, there is a good agreement between the crossflow model
((7) –(9)) and the experimental data for both the mixtures. On the other hand,
the perfect mixing model always underestimates the composition of the
retentate and the permeate streams leaving the membrane module. Therefore, the
optimization results obtained using (7)–(9) as the permeator model are more
reliable than those obtained using the perfect mixing model.

We are not aware of literature that provides, for liquid mixtures, detailed
experimental data regarding the composition of permeate and retentate streams
as a function of stage cut. We will nevertheless use (7)–(9) for liquid
mixtures because the agreement between the local flux determined using (3) and
the experimental value is very good [47, 48]. We recognize that this may not
sufficiently validate (7)–(9) for liquid mixtures, but unfortunately the
verification of composition profile is not possible in the absence of
experimental data.

Figure 3: Comparison between experimental data for the separation of O2/N2
[38] and the predictions obtained from crossflow model and the perfect mixing
model. Here, $P^{\textnormal{out}}/P^{\textnormal{per}}=7.83\text{
bar}/0.93\text{ bar}=8.4$, $x_{\ce{O2}}^{\textnormal{in}}=0.205$, membrane
permselectivity, $S=5.3$ (calculated based on the average values of component
permeances reported in [38]). Figure 4: Comparison between experimental data
for the separation of CO2/CH4 [49] and the predictions obtained from the
cross-flow model and the perfect mixing model. Here,
$P^{\textnormal{out}}/P^{\textnormal{per}}=4.05\text{ bar}/1.01\text{ bar}=4$,
$x_{\ce{CO2}}^{\textnormal{in}}=0.60$, the membrane permselectivity at 25 °C
is: 3.58, and the membrane permselectivity at 65 °C is: 2.9 [49].

###  3.2 Minimum Selectivity Requirement

By definition, selectivity is the ratio of the permeances of the more
permeable and less permeable components, so $S>1$. Now, consider the flux
equation in (5) in the following form

$y-x=k(S-1)^{2}\frac{y(1-y)}{S-(S-1)y}.$ (11)

Since $0\leq y\leq 1$, the sign of the RHS depends on the sign of $k$. For a
gaseous mixture, from (6), it can be shown that $k\geq 0$ regardless of the
value of the selectivity, provided $r\geq 1$. This implies $y-x\geq 0$, or the
local mole fraction on the permeate side is greater than that on the retentate
side. This is consistent with the expected behavior of the permeation process.
On the contrary, for a given liquid mixture (i.e., for a given $C_{A}$ and
$C_{B}$), $k$ can be negative even for some $S>1$. In that case, $y-x\leq 0$,
and thus, the permeate from the membrane module is enriched in the less
permeable component. This phenomenon, referred to as negative rejection [50,
51], has been observed experimentally in the context of separation of
phenol/water via reverse osmosis [50]. In this work, we avoid negative
rejection by choosing

$S>\frac{1-e^{-C_{B}\,u}}{1-e^{-C_{A}\,u}}.$ (12)

Then, from (6) and (5), $k$ and $y-x$ are guaranteed to be positive.

Further, in a typical permeation process, $y$ increases with increase in
$r$/$\Delta P^{trans}$ for a given $x$. Therefore, for a given $x$, we require

$\frac{dy}{dr}=\frac{dy}{du}\cdot\frac{du}{dr}\geq 0$ (13)

for a gaseous mixture and

$\frac{dy}{d(\Delta P^{trans})}=\frac{dy}{du}\cdot\frac{du}{d(\Delta
P^{trans})}\geq 0$ (14)

for a liquid mixture. Since $du/dr\geq 0$ and $du/d(\Delta P^{trans})\geq 0$
(see (2a)), we require $dy/du\geq 0$. Differentiating (5) with respect to $u$
yields

$\frac{dy}{du}=\frac{dk}{du}(S-1)^{2}\left[\frac{Sx}{y^{2}}+\frac{1-x}{(1-y)^{2}}\right]^{-1}.$
(15)

Since $0\leq y\leq 1$ and $0\leq x\leq 1$, $dy/du\geq 0$ only when $dk/du\geq
0$. For a gaseous mixture, from (6), it can be verified that $dk/du\geq 0$
regardless of the value of the selectivity. On the other hand, for a liquid
mixture, $dk/du\geq 0$ only when

$S>\frac{C_{B}}{C_{A}}e^{(C_{A}-C_{B})u}.$ (16)

We obtain the minimum selectivity needed by combining (12) and (16), and
noting that $u$ is a decision variable, as

$S>\max_{u^{\textnormal{lo}}\leq u\leq
u^{\textnormal{up}}}\left\\{\frac{1-e^{-C_{B}\,u}}{1-e^{-C_{A}\,u}},\frac{C_{B}}{C_{A}}e^{(C_{A}-C_{B})u}\right\\}.$
(17)

Here, $u^{\textnormal{lo}}=(\Delta P^{trans})^{\textnormal{lo}}$ and
$u^{\textnormal{up}}=(\Delta P^{trans})^{\textnormal{up}}$ (resp.
$u^{\textnormal{lo}}=\ln r^{\textnormal{lo}}$ and $u^{\textnormal{up}}=\ln
r^{\textnormal{up}}$), where $(\Delta P^{trans})^{\textnormal{lo}}$ and
$(\Delta P^{trans})^{\textnormal{up}}$ (resp. $r^{\textnormal{lo}}$ and
$r^{\textnormal{up}}$) denote the lower and upper bounds on the trans-membrane
pressure difference (resp. pressure ratio) for the given liquid (resp.
gaseous) mixture. Note that, (17) reduces to $S>1$ for a gaseous mixture. On
the other hand, the minimum selectivity needed to separate the given liquid
mixture depends on the molar volumes of the constituent components ($C_{A}$
and $C_{B}$) and on the admissible range of the trans-membrane pressure
difference. In this work, we consider only those cases where (17) holds.

###  3.3 Properties of the Permeator Model

Here, we describe the properties of the permeator model which we will use in
§4.5 to derive additional cuts to the MINLP. Provided $\eqref{eq:min-Selec}$
holds,

  * P1

The mole fraction of component $A$ on the retentate side decreases along the
length of the membrane module. Mathematically, this can be readily shown using
(4). Since $y-x\geq 0$ (see §3.2) and $f\geq 0$, $dx/df\geq 0$. Further, since
$f$ decreases along the length of the membrane module, $x$ also decreases.
This implies that $x^{\textnormal{out}}\leq x^{\textnormal{in}}$.

  * P2

The mole fraction of component $A$ in the permeate is at least as high as the
mole fraction in the feed. Mathematically, from (10), $x^{\textnormal{in}}$
can be expressed as a convex combination of $x^{\textnormal{out}}$ and
$y^{\textnormal{per}}$ i.e.,
$x^{\textnormal{in}}=(1-\theta)x^{\textnormal{out}}+\theta
y^{\textnormal{per}}$. Since $x^{\textnormal{out}}\leq x^{\textnormal{in}}$
from P1, $y^{\textnormal{per}}\geq x^{\textnormal{in}}$.

  * P3

The local mole fraction of component $A$ on the permeate side ($y$) increases
monotonically with an increase in the local mole fraction on the retentate
side ($x$). This is because

$\displaystyle\frac{dy}{dx}=y\left[x+\frac{k(S-1)^{2}y^{2}}{[S-(S-1)y]^{2}}\right]^{-1}$
(18)

(see §A.3 for derivation) is non-negative.

  * P4

The local mole fraction of component $A$ on the permeate side ($y$) decreases
along the length of the membrane module. This is inferred from properties P1
and P3. Thus, $y^{\textnormal{out}}\leq y^{\textnormal{in}}$. Further, since
the overall permeate is the aggregate of permeate from differential membrane
elements, its composition $(y^{\textnormal{per}})$ satisfies
$y^{\textnormal{out}}\leq y^{\textnormal{per}}\leq y^{\textnormal{in}}$.

##  4 Problem Formulation

###  4.1 Problem statement

The required input to the problem consists of (i) the molar flowrate and the
composition of the feed and the product streams, (ii) efficiencies of
compressors and turbines (for gas separations)/pumps and turbocharger (for
liquid separations), (iii) membrane perm-selectivity, (iv) the range of
admissible operating pressure ratio/trans-membrane pressure difference, (v)
temperature of the mixture, and (vi) molar volume of the constituent
components for liquid mixtures (see Table 1).

Given a binary mixture along with all the required inputs, the problem is then
to identify the membrane cascade that requires at most $N$ stages and consumes
least power for the separation.

Table 1: List of input parameters Symbol |  Definition
---|---
$N$ |  Maximum number of stages in the cascade
$F$, $F^{\textnormal{per}}$, $F^{\textnormal{out}}$ |  Molar flowrate of the given feed mixture, permeate product, and retentate product streams, respectively
$X^{F}$, $Y^{\textnormal{per}}$, $X^{\textnormal{out}}$ |  Mole fraction of component $A$ in the feed mixture, permeate product, and retentate product streams, respectively
$\eta^{comp}$ |  Isothermal compressor efficiency
$\eta^{pump}$ |  Pump efficiency
$\eta^{TC}$ |  Turbocharger efficiency
$V_{A},V_{B}$ |  Liquid molar volume of components $A$ and $B$, respectively (needed only for liquid mixtures)
$V^{F}$ |  Liquid molar volume of the feed mixture, calculated as $V^{F}=X^{F}V_{A}+(1-X^{F})V_{B}$ (needed only for liquid mixtures)
$V^{\textnormal{out}}$ |  Liquid molar volume of the retentate product, calculated as $V^{\textnormal{out}}=X^{\textnormal{out}}V_{A}+(1-X^{\textnormal{out}})V_{B}$ (needed only for liquid mixtures)
$T$ |  Absolute temperature of the feed mixture
$S$ |  Membrane permselectivity ($\mathit{PM}_{A}/\mathit{PM}_{B}$ )
$[r^{\textnormal{lo}},r^{\textnormal{up}}]$ |  Admissible range of trans-membrane pressure ratio for a gaseous mixture
$[(\Delta P^{trans})^{\textnormal{lo}},(\Delta P^{trans})^{\textnormal{up}}]$ |  Admissible range of trans-membrane pressure difference for a liquid mixture

###  4.2 Membrane Cascade Superstructure

Figure 5(a) (resp. Figure 5(b)) shows the superstructure which embeds cascades
requiring at most $N$ stages for the separation of a gaseous mixture (resp.
liquid mixture). The splitter $\mathcal{F}$ (see Figure 5) splits the feed
mixture into $N$ streams which are sent to mixers $\mathcal{M}_{1}$ through
$\mathcal{M}_{N}$. Each mixer $\mathcal{M}_{j}$, $j=1,\dots,N-1$, supplies the
feed to stage $j$ after mixing the retentate from the stage $j-1$ and the
streams from the splitters $\mathcal{F}$, $\mathcal{S}_{j+1}$ and
$\mathcal{S}_{j+2}$ (see Figure 5). Mixer $\mathcal{M}_{N}$ supplies the feed
to stage $N$ after mixing the streams from the splitters $\mathcal{F}$ and
$\mathcal{O}_{N-1}$. The permeate from each stage $j$ is sent to the splitter
$\mathcal{S}_{j}$. Each splitter $S_{j}$, $j=3,\dots,N$, splits the stream
into two streams which are sent to mixers $\mathcal{M}_{j-1}$ and
$\mathcal{M}_{j-2}$. Splitter $\mathcal{S}_{2}$ splits the permeate from stage
2 into two streams which are sent to mixers $\mathcal{M}_{1}$ and
$\mathcal{P}$. Splitter $\mathcal{S}_{1}$ sends the permeate from stage 1 to
the mixer $\mathcal{P}$. On the other hand, the retentate from each stage $j$,
$j=1,\dots,N-2$, is sent to the mixer $\mathcal{M}_{j}$. The retentate from
the stage $N-1$ (resp. $N$) is sent to the splitter $\mathcal{O}_{N-1}$ (resp.
$\mathcal{O}_{N}$). Splitter $\mathcal{O}_{N-1}$ splits the inlet stream into
two streams which are sent to the mixers $\mathcal{M}_{N}$ and $\mathcal{R}$,
and the splitter $\mathcal{O}_{N}$ sends the retentate from the stage $N$ to
the mixer $\mathcal{R}$. Mixer $\mathcal{P}$ (resp. $\mathcal{R}$) mixes the
streams from the splitters $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ (resp.
$\mathcal{O}_{N-1}$ and $\mathcal{O}_{N}$) and produces the permeate product
stream (resp. retentate product stream) (see Figure 5).

Figure 5: Superstructure for the separation of a) gaseous mixtures b) liquid
mixtures

In this work, we assume that the given gaseous mixture is at a high pressure
and the products are also desired at a high pressure. However, if the gaseous
mixture is not available at a high pressure, then an additional compressor can
be included in the superstructure to compress the feed. Further, if the
products are not desired at a high pressure, then (i) the compressor used for
compressing the permeate product stream (see Figure 5(a)) can be eliminated
from the superstructure, and (ii) a turbine can be included in the
superstructure to recover work by expanding the retentate product stream. On
the other hand, we assume that the given liquid mixture is always available at
a low pressure. We increase the pressure of the feed mixture to an
intermediate value using a pump. Next, we send the feed mixture to a
turbocharger where it is further pressurized to the desired pressure by
transferring the work recovered from the expansion of the retentate product
stream (see Figure 5(b)).

Further, we impose the following restrictions.

  1. 1.

All the splitters
($\mathcal{F},\mathcal{S}_{1},\dots,\mathcal{S}_{N},\mathcal{O}_{N-1},$ and
$\mathcal{O}_{N}$) direct the material flow to one arc entirely. In other
words, only one of the arcs arising from a splitter can contain nonzero
material flow.

  2. 2.

The trans-membrane pressure ratio/trans-membrane pressure difference is the
same for all the stages in the cascade.

Under the above restrictions, we believe that the optimal solution satisfies
the following property.

  * P5

In the optimal solution, the mole fraction of the most permeable component in
the inlet, retentate and permeate streams decrease from stage 1 through stage
$N$ i.e., $x_{1}^{\textnormal{in}}\geq\dots\geq x_{N}^{\textnormal{in}}$,
$x_{1}^{\textnormal{out}}\geq\dots\geq x_{N}^{\textnormal{out}}$, and
$y_{1}^{\textnormal{per}}\geq\dots\geq y_{N}^{\textnormal{per}}$.

The above property is based on physical intuition and empirical observation.
We use this property in §4.5 to derive additional cuts to the MINLP.

###  4.3 Objective function

Figure 6 shows the variables used in our formulation and their definition can
be found in Table 2.

Figure 6: Variables used in the optimization formulation. a) single membrane module, b) splitter $\mathcal{S}_{2}$, c) splitter $\mathcal{O}_{N-1}$ Table 2: List of variables Symbol |  Definition
---|---
$u$ |  Auxiliary variable defined in (2a)
$k$ |  Auxiliary variable defined in (W16)
$f^{F}_{j}$ |  Molar flow rate along the arc connecting $\mathcal{F}$ and $\mathcal{M}_{j}$
$f^{\textnormal{in}}_{j}$ |  Molar flow rate entering stage $j$ from $\mathcal{M}_{j}$
$f^{\textnormal{out}}_{j}$ |  Molar flow rate of the retentate stream leaving stage $j$
$f^{\textnormal{out},1}_{j}$ |  Molar flow rate along the arc connecting $\mathcal{M}_{j+1}$ and $\mathcal{O}_{j}$
$f^{\textnormal{out,bypass}}_{j}$ |  Molar flow rate along the arc connecting $\mathcal{R}$ and $\mathcal{O}_{j}$
$f^{\textnormal{per}}_{j}$ |  Molar flow rate of the permeate stream leaving stage $j$
$f^{\textnormal{recy},p}_{j}$ |  Molar flow rate along the arc connecting $\mathcal{M}_{j-p}$ and $\mathcal{S}_{j}$, $p=1,2$
$f^{\textnormal{per,bypass}}_{j}$ |  Molar flow rate along the arc connecting $\mathcal{P}$ and $\mathcal{S}_{j}$
$x_{j}^{\textnormal{in}}$ |  Mole fraction of component $A$ in the inlet stream to stage $j$
$x_{j}^{\textnormal{out}}$ |  Mole fraction of component $A$ in the retentate stream leaving stage $j$
$y^{\textnormal{per}}_{j}$ |  Mole fraction of component $A$ in the permeate stream leaving stage $j$
$y^{\textnormal{in}}_{j}$ |  Local mole fraction of component $A$ at the entrance of stage $j$ on permeate side
$y^{\textnormal{out}}_{j}$ |  Local mole fraction of component $A$ at the exit of stage $j$ on permeate side
$z_{j}^{\textnormal{in}},\;z_{j}^{\textnormal{out}}$ |  Auxiliary variables defined in (W18)
$\theta_{j}$ |  Stage cut of stage $j$
$D_{j}$ |  Auxiliary variable defined in (23)
$\omega^{F}_{j}$ |  Binary variable introduced to regulate the flow along the arc connecting $\mathcal{F}$ and $\mathcal{M}_{j}$
$\omega^{\textnormal{out},1}_{j}$ |  Binary variable introduced to regulate the flow along the arc connecting $\mathcal{M}_{j+1}$ and $\mathcal{O}_{j}$
$\omega^{\textnormal{out,bypass}}_{j}$ |  Binary variable introduced to regulate the flow along the arc connecting $\mathcal{R}$ and $\mathcal{O}_{j}$
$\omega^{\textnormal{recy},p}_{j}$ |  Binary variable introduced to regulate the flow along the arc connecting $\mathcal{M}_{j-p}$ and $\mathcal{S}_{j}$, $p=1,2$
$\omega^{\textnormal{per,bypass}}_{j}$ |  Binary variable introduced to regulate the flow along the arc connecting $\mathcal{P}$ and $\mathcal{S}_{j}$

First, consider the separation of a gaseous mixture. We assume that the given
gas mixture is at a high pressure, and that the products are needed at the
same pressure as the feed. In this case, energy is needed to compress the
outlet streams from splitters $\mathcal{S}_{j}$ to $r-$times its pressure (see
Figure 5(a)). We estimate the compressor power using an isothermal compressor
efficiency, $\eta^{comp}$. The net power required for the cascade is obtained
as

$W_{G}=\frac{R\,T}{\eta^{comp}}\left(\sum_{j=2}^{N}f_{j}^{\textnormal{recy},1}+\sum_{j=3}^{N}f_{j}^{\textnormal{recy},2}+F^{\textnormal{per}}\right)\ln
r.$ (19)

Next, consider the separation of a liquid mixture. We assume that the feed
mixture and the product streams are at pressure $P^{\textnormal{per}}$. In
this case, energy is needed to pump the feed mixture from
$P^{\textnormal{per}}$ to an intermediate pressure $P^{TC}$, and to pump the
outlet streams from each splitter $\mathcal{S}_{j}$, $j=2,\dots,N$, from
$P^{\textnormal{per}}$ to $P^{\textnormal{out}}$ (see Figure 5(b)). Assuming a
pump efficiency of $\eta^{pump}$, the net power required for the cascade is
obtained as

$\displaystyle
W_{L}=\underbrace{\frac{V^{F}}{\eta^{pump}}F\,(P^{TC}-P^{\textnormal{per}})}_{\text{to
pump the
feed}}+\underbrace{\frac{1}{\eta^{pump}}\left[\sum_{j=2}^{N}f_{j}^{\textnormal{recy},1}V^{\textnormal{per}}_{j}\,(P^{\textnormal{out}}-P^{\textnormal{per}})+\sum_{j=3}^{N}f_{j}^{\textnormal{recy},2}V^{\textnormal{per}}_{j}(P^{\textnormal{out}}-P^{\textnormal{per}})\right]}_{\text{to
pump the recycle streams}}.$ (20)

Here, $V^{F}$ and $V_{j}^{\textnormal{per}}$ denote the molar volume of the
feed and the molar volume of the permeate stream from stage $j$, respectively.
As mentioned earlier, the turbocharger pressurizes the feed from $P^{TC}$ to
$P^{\textnormal{out}}$ while expanding the retentate product stream from
$P^{\textnormal{out}}$ to $P^{\textnormal{per}}$. Let, $\eta^{TC}$ denote the
efficiency of the turbocharger. Then,

$\displaystyle
V^{F}\,F\,(P^{\textnormal{out}}-P^{TC})=\eta^{TC}\,\left[V^{\textnormal{out}}\,F^{\textnormal{out}}\,(P^{\textnormal{out}}-P^{\textnormal{per}})\right].$
(21)

Here, the LHS corresponds to the power needed to pressurize the feed mixture,
the term inside the square brackets on the RHS corresponds to the power
produced from the expansion of the retentate product stream, and
$V^{\textnormal{out}}$ denotes the molar volume of the retentate product
stream. We substitute $P^{TC}$ from (21) in (20) to obtain

$\displaystyle W_{L}=\frac{1}{\eta^{pump}}\left(FV^{F}\,\Delta
P^{trans}-\eta^{TC}F^{\textnormal{out}}V^{\textnormal{out}}\,\Delta
P^{trans}+\sum_{j=2}^{N}f_{j}^{\textnormal{recy},1}V^{\textnormal{per}}_{j}\,\Delta
P^{trans}+\sum_{j=3}^{N}f_{j}^{\textnormal{recy},2}V^{\textnormal{per}}_{j}\,\Delta
P^{trans}\right).$ (22)

From (19) and (22), we obtain a unified objective function that is applicable
for both gaseous and liquid mixtures using (2a) as
$D_{0}u+D_{1}F^{\textnormal{per}}\cdot
u+\sum_{j=2}^{N}D_{j}f_{j}^{\textnormal{recy},1}\cdot
u+\sum_{j=3}^{N}D_{j}f_{j}^{\textnormal{recy},2}\cdot
u+D_{N+1}F^{\textnormal{out}}\cdot u$, where

$\displaystyle D_{0}=\begin{cases}0,&\text{for a gaseous mixture}\\\
\frac{V^{F}\,F}{\eta^{pump}},&\text{for a liquid mixture}\end{cases}$ (23a)
$\displaystyle D_{1}=\begin{cases}\frac{R\,T}{\eta^{comp}},&\text{for a
gaseous mixture}\\\ 0,&\text{for a liquid mixture}\end{cases}$ (23b)
$\displaystyle D_{j}=\begin{cases}\frac{R\,T}{\eta^{comp}},&\text{for a
gaseous mixture}\\\ \frac{V_{j}^{\textnormal{per}}}{\eta^{pump}},&\text{for a
liquid mixture}\end{cases},\quad j=2,\dots,N$ (23c) $\displaystyle
D_{N+1}=\begin{cases}0,&\text{for a gaseous mixture}\\\
-\frac{V^{\textnormal{out}}\,\eta^{TC}}{\eta^{pump}},&\text{for a liquid
mixture}\end{cases}$ (23d)

By assuming ideal mixing, we obtain $V^{F}=V_{A}X^{F}+V_{B}(1-X^{F})$,
$V_{j}^{\textnormal{per}}=V_{A}y_{j}^{\textnormal{per}}+V_{B}(1-y_{j}^{\textnormal{per}})$,
and
$V^{\textnormal{out}}=V_{A}X^{\textnormal{out}}+V_{B}(1-X^{\textnormal{out}})$,
where $V_{A}$ and $V_{B}$ denote the molar volume of the pure components $A$
and $B$, respectively. Note that $D_{0},\dots,D_{N+1}$ are parameters (resp.
variables) for a gaseous (resp. liquid) mixture.

###  4.4 MINLP formulation

Here, we present our mixed-integer nonlinear program (W) for identifying the
optimal membrane cascade requiring at most $N$ membrane stages. Consider
Figure 6. Let $\mathcal{J}=\\{1,\dots,N\\}$. In the following, acronyms OMB
and CMB stand for overall mass balance and mass balance of component $A$,
$\delta_{(\cdot)}=\\{1,\text{ if }(\cdot)\text{ is true};\;0,\text{
otherwise}\\}$, $(\cdot)^{\textnormal{lo}}$ and $(\cdot)^{\textnormal{up}}$
denote the lower and upper bounds on $(\cdot)$.

$\displaystyle(\text{W}):\;\;\min\quad$ $\displaystyle
D_{0}u+D_{1}F^{\textnormal{per}}\cdot
u+\sum_{j=2}^{N}D_{j}f_{j}^{\textnormal{recy},1}\cdot u$
$\displaystyle+\sum_{j=3}^{N}D_{j}f_{j}^{\textnormal{recy},2}\cdot
u+D_{N+1}F^{\textnormal{out}}\cdot u,$ (Objective Function) (W1)
$\displaystyle\text{s.t.},\quad$ $\displaystyle\eqref{eq:Dvar-defs},$
$\displaystyle(\text{Definition of $D_{j}$})$ (W2)
$\displaystyle\sum_{j=1}^{N}f_{j}^{F}=F,$ $\displaystyle(\text{OMB around
}\mathcal{F})$ (W3) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle\left.\begin{aligned}
f_{j}^{\textnormal{in}}=&f_{j-1}^{\textnormal{out}}\delta_{2\leq j\leq
N-1}+f_{j-1}^{\textnormal{out},1}\delta_{j=N}+f_{j+1}^{\textnormal{recy},1}\delta_{j\leq
N-1}\\\ &+f_{j+2}^{\textnormal{recy},2}\delta_{j\leq
N-2}+f_{j}^{F}\end{aligned}\right\\},$ $\displaystyle(\text{OMB around
}\mathcal{M}_{j})$ (W4) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle\left.\begin{aligned}
&f_{j}^{\textnormal{in}}x_{j}^{\textnormal{in}}=f_{j-1}^{\textnormal{out}}x_{j-1}^{\textnormal{out}}\delta_{2\leq
j\leq
N-1}+f_{j-1}^{\textnormal{out},1}x_{j-1}^{\textnormal{out}}\delta_{j=N}\\\
&+f_{j+1}^{\textnormal{recy},1}y^{\textnormal{per}}_{j+1}\delta_{j\leq
N-1}+f_{j+2}^{\textnormal{recy},2}y^{\textnormal{per}}_{j+2}\delta_{j\leq
N-2}+f_{j}^{F}X^{F}\end{aligned}\right\\},$ $\displaystyle(\text{CMB around
}\mathcal{M}_{j})$ (W5) $\displaystyle\left.\begin{aligned}
&f^{\textnormal{out}}_{N-1}=f_{N-1}^{\textnormal{out},1}+f_{N-1}^{\textnormal{out,bypass}}\\\
&f^{\textnormal{out}}_{N}=f_{N}^{\textnormal{out,bypass}}\end{aligned}\right\\},$
$\displaystyle(\text{OMB around }\mathcal{O}_{j})$ (W6)
$\displaystyle\left.\begin{aligned}
&f^{\textnormal{out}}_{N-1}x^{\textnormal{out}}_{N-1}=f_{N-1}^{\textnormal{out},1}x^{\textnormal{out}}_{N-1}+f_{N-1}^{\textnormal{out,bypass}}x^{\textnormal{out}}_{N-1}\\\
&f^{\textnormal{out}}_{N}x^{\textnormal{out}}_{N}=f_{N}^{\textnormal{out,bypass}}x^{\textnormal{out}}_{N}\end{aligned}\right\\},$
$\displaystyle(\text{CMB around }\mathcal{O}_{j})$ (W7)
$\displaystyle\forall\;j\in\mathcal{J},\quad$ $\displaystyle
f^{\textnormal{per}}_{j}=f_{j}^{\textnormal{recy},1}\delta_{j\geq
2}+f_{j}^{\textnormal{recy},2}\delta_{j\geq
3}+f_{j}^{\textnormal{per,bypass}}\delta_{j\leq 2},$ $\displaystyle(\text{OMB
around }\mathcal{S}_{j})$ (W8) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle\left.\begin{aligned}
f^{\textnormal{per}}_{j}y^{\textnormal{per}}_{j}=&f_{j}^{\textnormal{recy},1}y^{\textnormal{per}}_{j}\delta_{j\geq
2}+f_{j}^{\textnormal{recy},2}y^{\textnormal{per}}_{j}\delta_{j\geq 3}\\\
&+f_{j}^{\textnormal{per,bypass}}y^{\textnormal{per}}_{j}\delta_{j\leq
2}\end{aligned}\right\\},$ $\displaystyle(\text{CMB around }\mathcal{S}_{j})$
(W9) $\displaystyle\forall\;j\in\mathcal{J},\quad$ $\displaystyle
f_{j}^{\textnormal{in}}=f_{j}^{\textnormal{per}}+f_{j}^{\textnormal{out}},$
$\displaystyle(\text{OMB around stage }j)$ (W10)
$\displaystyle\forall\;j\in\mathcal{J},\quad$ $\displaystyle
f_{j}^{\textnormal{in}}x_{j}^{\textnormal{in}}=f_{j}^{\textnormal{per}}y_{j}^{\textnormal{per}}+f_{j}^{\textnormal{out}}x_{j}^{\textnormal{out}},$
$\displaystyle(\text{CMB around stage }j)$ (W11)
$\displaystyle\sum_{j=N-1}^{N}f_{j}^{\textnormal{out,bypass}}=F^{\textnormal{out}}$
$\displaystyle(\text{OMB around }\mathcal{R})$ (W12)
$\displaystyle\sum_{j=N-1}^{N}f_{j}^{\textnormal{out,bypass}}x^{\textnormal{out}}_{j}=F^{\textnormal{out}}X^{\textnormal{out}}$
$\displaystyle(\text{CMB around }\mathcal{R})$ (W13)
$\displaystyle\sum_{j=1}^{2}f_{j}^{\textnormal{per,bypass}}=F^{\textnormal{per}}$
$\displaystyle(\text{OMB around }\mathcal{P})$ (W14)
$\displaystyle\sum_{j=1}^{2}f_{j}^{\textnormal{per,bypass}}y^{\textnormal{per}}_{j}=F^{\textnormal{per}}Y^{\textnormal{per}}$
$\displaystyle(\text{CMB around }\mathcal{P})$ (W15)
$\displaystyle(S-1)^{2}k=(S-1)-\left(Se^{-C_{A}u}-e^{-C_{B}u}\right),$
$\displaystyle(\text{Definition of }k)$ (W16)
$\displaystyle\forall\;j\in\mathcal{J},\quad$ $\displaystyle
f_{j}^{\textnormal{per}}=\theta_{j}f_{j}^{\textnormal{in}},$ (Definition of
stage cut) (W17) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle\left.\begin{aligned} &S\ln y_{j}^{\textnormal{out}}-S\ln
y_{j}^{\textnormal{in}}-\ln(1-y_{j}^{\textnormal{out}})+\ln(1-y_{j}^{\textnormal{in}})\\\
&-(S-1)^{2}(k\ln z_{j}^{\textnormal{out}})+(S-1)^{2}(k\ln
z_{j}^{\textnormal{in}})\\\ &=(S-1)^{2}\left[k\ln(1-\theta_{j})\right],\\\
&z_{j}^{\textnormal{in}}=y_{j}^{\textnormal{in}}-x_{j}^{\textnormal{in}},\\\
&z_{j}^{\textnormal{out}}=y_{j}^{\textnormal{out}}-x_{j}^{\textnormal{out}},\\\
&k+(S-1)(ky_{j}^{\textnormal{in}})-S\left[\frac{k}{S-(S-1)y_{j}^{\textnormal{in}}}\right]=z_{j}^{\textnormal{in}},\\\
&k+(S-1)(ky_{j}^{\textnormal{out}})-S\left[\frac{k}{S-(S-1)y_{j}^{\textnormal{out}}}\right]=z_{j}^{\textnormal{out}},\end{aligned}\right\\}$
$\displaystyle(\text{Permeator Model})$ (W18)
$\displaystyle\left.\begin{aligned}
&\omega_{j}^{F}&&\in\\{0,1\\},&&\forall\;j\in\\{1,\dots,N\\}\\\
&\omega_{j}^{\textnormal{recy},1}&&\in\\{0,1\\},&&\forall\;j\in\\{2,\dots,N\\}\\\
&\omega_{j}^{\textnormal{recy},2}&&\in\\{0,1\\},&&\forall\;j\in\\{3,\dots,N\\}\\\
&\omega_{j}^{\textnormal{per,bypass}}&&\in\\{0,1\\},&&\forall\;j\in\\{1,2\\}\\\
&\omega_{N-1}^{\textnormal{out,bypass}}&&\in\\{0,1\\}\\\
&\omega_{N-1}^{\textnormal{out},1}&&\in\\{0,1\\}\end{aligned}\right\\}$ (W19)
$\displaystyle\sum_{j=1}^{N}\omega_{j}^{F}=1,$ (W20)
$\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle\omega^{\textnormal{recy},1}_{j}\delta_{j\geq
2}+\omega^{\textnormal{recy},2}_{j}\delta_{j\geq
3}+\omega^{\textnormal{per,bypass}}_{j}\delta_{j\leq 2}=1$ (W21)
$\displaystyle\left.\begin{aligned}
&\omega_{N-1}^{\textnormal{out},1}+\omega_{N-1}^{\textnormal{out,bypass}}=1,\\\
&\omega_{N}^{\textnormal{out,bypass}}=1\end{aligned}\right\\},$ (W22)
$\displaystyle\left.\begin{aligned} &0\leq
f_{j}^{F}\leq\omega_{j}^{F}(f_{j}^{F})^{\textnormal{up}},&&\forall\;j\in\\{1,\dots,N\\}\\\
&0\leq
f_{j}^{\textnormal{recy},1}\leq\omega_{j}^{\textnormal{recy},1}(f_{j}^{\textnormal{recy},1})^{\textnormal{up}},&&\forall\;j\in\\{2,\dots,N\\}\\\
&0\leq
f_{j}^{\textnormal{recy},2}\leq\omega_{j}^{\textnormal{recy},2}(f_{j}^{\textnormal{recy},2})^{\textnormal{up}},&&\forall\;j\in\\{3,\dots,N\\}\\\
&0\leq
f_{j}^{\textnormal{per,bypass}}\leq\omega_{j}^{\textnormal{per,bypass}}(f_{j}^{\textnormal{per,bypass}})^{\textnormal{up}},&&\forall\;j\in\\{1,2\\}\\\
&0\leq
f_{j}^{\textnormal{out,bypass}}\leq\omega_{j}^{\textnormal{out,bypass}}(f_{j}^{\textnormal{out,bypass}})^{\textnormal{up}},&&\forall\;j\in\\{N-1,N\\}\\\
&0\leq
f_{N-1}^{\textnormal{out},1}\leq\omega_{N-1}^{\textnormal{out},1}(f_{N-1}^{\textnormal{out},1})^{\textnormal{up}}\end{aligned}\right\\}$
(W23) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle\left.\begin{aligned}
&0\leq(\cdot)\leq(\cdot)^{\textnormal{up}},\quad\forall\;(\cdot)\in\left\\{f_{j}^{\textnormal{in}},f_{j}^{\textnormal{out}},f_{j}^{\textnormal{per}}\right\\}\\\
&(\cdot)^{\textnormal{lo}}\leq(\cdot)\leq(\cdot)^{\textnormal{up}},\quad\forall\;(\cdot)\in\left\\{u,k,\theta_{j},x_{j}^{\textnormal{in}},x_{j}^{\textnormal{out}},\right.\\\
&\hskip
136.57323pt\left.y_{j}^{\textnormal{per}},y_{j}^{\textnormal{in}},y_{j}^{\textnormal{out}},z_{j}^{\textnormal{in}},z_{j}^{\textnormal{out}}\right\\}\end{aligned}\right\\}$
(Bounds on variables) (W24)

The formulation of the objective function is described in the previous
subsection. We now describe the formulation of constraints.

Mass balance constraints: (W3) models overall mass balance around the feed
splitter $\mathcal{F}$ (see the superstructures in Figure 5). (W4), (W6),
(W8), and (W10) (resp. (W5), (W7), (W9), and (W11)) model overall mass balance
(resp. mass balance of component $A$) around mixer $\mathcal{M}_{j}$,
splitters $\mathcal{O}_{j}$ and $\mathcal{S}_{j}$, and membrane stage $j$,
respectively. Mass balances of component $B$ are implied from the overall mass
balances and the mass balances on component $A$, so we do not impose them
explicitly. However, although (W7) (resp. (W9)) is implied from (W6) (resp.
(W8)), we impose it explicitly, because it is not implied in the relaxation
where the bilinear terms appear in relaxed form. Next, (W12) and (W13) (resp.
(W14) and (W15)) model the overall mass balance and the mass balance of
component $A$ around the mixer $\mathcal{R}$ (resp. $\mathcal{P}$).

Permeator model constraints: (W16) is the same as (6) and it computes the
value of $k$. (W17) computes the value of the stage cut. (W18) is the
permeator model in (7)–(9). Observe that we introduced auxiliary variables
$z_{j}^{\textnormal{in}}$ and $z_{j}^{\textnormal{out}}$ for
$y_{j}^{\textnormal{in}}-x_{j}^{\textnormal{in}}$ and
$y_{j}^{\textnormal{out}}-x_{j}^{\textnormal{out}}$ along with the constraints
$(z_{j}^{\textnormal{in}})^{\textnormal{lo}}\leq z_{j}^{\textnormal{in}}$ and
$(z_{j}^{\textnormal{out}})^{\textnormal{lo}}\leq z_{j}^{\textnormal{out}}$
(see (W24)), where $0<(z_{j}^{\textnormal{in}})^{\textnormal{lo}}$ and
$0<(z_{j}^{\textnormal{out}})^{\textnormal{lo}}$. The choice of the lower
bounds will be discussed shortly. Without the auxiliary variables and the
bound constraints, BARON reports an error since it cannot infer
$y_{j}^{\textnormal{in}}-x_{j}^{\textnormal{in}}>0$ and
$y_{j}^{\textnormal{out}}-x_{j}^{\textnormal{out}}>0$, which are needed to
well-define $\ln(y_{j}^{\textnormal{in}}-x_{j}^{\textnormal{in}})$ and
$\ln(y_{j}^{\textnormal{out}}-x_{j}^{\textnormal{out}})$ terms (see (7)).
Further, we have disaggregated all log terms i.e., each log term of the form
$\ln(y_{j}^{\textnormal{out}}/y_{j}^{\textnormal{in}})$ is expressed as $\ln
y_{j}^{\textnormal{out}}-\ln y_{j}^{\textnormal{in}}$. Without disaggregation,
a typical factorable relaxation procedure first introduces an auxiliary
variable for the fraction $y_{j}^{\textnormal{out}}/y_{j}^{\textnormal{in}}$,
and then relaxes the log term over the range of
$y_{j}^{\textnormal{out}}/y_{j}^{\textnormal{in}}$. Empirically, we observed
that BARON either fails to solve the MINLP, or does so slowly without the
disaggregation of the log terms.

Restricting flows along specific arcs: As mentioned in §4.2, in this work, we
require that all the splitters
($\mathcal{F},\mathcal{S}_{1},\dots,\mathcal{S}_{N},\mathcal{O}_{N-1},$ and
$\mathcal{O}_{N}$) direct the material flow to one arc entirely. (W19)–(W23)
model this requirement. First, we define binary variables (see (W19) for the
domain of index $j$) (i) $\omega_{j}^{F}=$ {1, if the splitter $\mathcal{F}$
directs the material flow to mixer $\mathcal{M}_{j}$; 0, otherwise}, (ii)
$\omega_{j}^{\textnormal{recy},1}=$ {1, if the splitter $\mathcal{S}_{j}$
directs the material flow to mixer $\mathcal{M}_{j-1}$; 0, otherwise}, (iii)
$\omega_{j}^{\textnormal{recy},2}=$ {1, if the splitter $\mathcal{S}_{j}$
directs the material flow to mixer $\mathcal{M}_{j-2}$; 0, otherwise}, (iv)
$\omega_{j}^{\textnormal{per,bypass}}=$ {1, if the spitter $\mathcal{S}_{j}$
directs the material flow to mixer $\mathcal{P}$; 0, otherwise}, (v)
$\omega_{j}^{\textnormal{out,bypass}}=$ {1, if the spitter $\mathcal{O}_{j}$
directs the material flow to mixer $\mathcal{R}$; 0, otherwise}, and (vi)
$\omega_{N-1}^{\textnormal{out},1}=$ {1, if the spitter $\mathcal{O}_{N-1}$
directs the material flow to mixer $\mathcal{M}_{N}$; 0, otherwise}. Next,
(W20)–(W22) ensure that only one arc, among all the arcs that are leaving from
a splitter, is chosen for directing the material flow. (W23) suppresses
material flow along the arcs when the corresponding binary variable takes the
value zero. Here, we choose $(f_{j}^{F})^{\textnormal{up}}=F$ because the
total flowrate along the arcs connecting splitter $\mathcal{F}$ and mixer
$\mathcal{M}_{j}$ cannot exceed the flowrate of the feed. Similarly, since the
total flowrate along the arcs connecting $\mathcal{S}_{j}$ and $\mathcal{P}$
(resp. $\mathcal{O}_{j}$ and $\mathcal{R}$) cannot exceed the flowrate of the
permeate product (resp. retentate product), we choose
$(f_{j}^{\textnormal{per,bypass}})^{\textnormal{up}}=F^{\textnormal{per}}$
(resp.
$(f_{j}^{\textnormal{out,bypass}})^{\textnormal{up}}=F^{\textnormal{out}}$).
However, a natural upper bound does not exist on the remaining flow variables.
Therefore, we choose a sufficiently large number for
$(f_{j}^{\textnormal{recy},1})^{\textnormal{up}},(f_{j}^{\textnormal{recy},2})^{\textnormal{up}},(f_{N-1}^{\textnormal{out},1})^{\textnormal{up}}$.

Bounds on variables: It is essential to have finite bounds on all variables,
especially those that are involved in nonlinear terms, in order to construct a
valid convex relaxation. (W23) bounds a few flow variables. (W24) bounds the
remaining variables in the problem. As before, there is no natural upper bound
on $f_{j}^{\textnormal{in}}$, $f_{j}^{\textnormal{per}}$ and
$f_{j}^{\textnormal{out}}$ variables, so we choose a sufficiently large number
for $(f_{j}^{\textnormal{in}})^{\textnormal{up}}$,
$(f_{j}^{\textnormal{per}})^{\textnormal{up}}$ and
$(f_{j}^{\textnormal{out}})^{\textnormal{up}}$.

For a gaseous (resp. liquid) mixture, let the admissible range of operating
pressure ratio (resp. trans-membrane pressure difference) for the chosen
membrane be $[r^{\textnormal{lo}},r^{\textnormal{up}}]$ (resp. [$(\Delta
P^{trans})^{\textnormal{lo}},(\Delta P^{trans})^{\textnormal{up}}$]). Then, we
choose the following as the lower and upper bounds on $u$.

$\displaystyle u^{\textnormal{lo}}=\begin{cases}\ln
r^{\textnormal{lo}},&\text{for a gaseous mixture}\\\ (\Delta
P^{trans})^{\textnormal{lo}},&\text{for a liquid mixture}\end{cases}$ (24a)
$\displaystyle u^{\textnormal{up}}=\begin{cases}\ln
r^{\textnormal{up}},&\text{for a gaseous mixture}\\\ (\Delta
P^{trans})^{\textnormal{up}},&\text{for a liquid mixture}\end{cases}$ (24b)

We obtain the lower and upper bounds on $k$ by substituting
$u=u^{\textnormal{lo}}$ and $u=u^{\textnormal{up}}$ in (6), respectively.
Since the chosen perm-selectivity is such that $k$ is a monotonically
increasing function of $u$ over the interval
$[u^{\textnormal{lo}},u^{\textnormal{up}}]$ (see §3.2), the choice of bounds
on $k$ is justified.

By definition, stage cut is the fraction of the total feed permeating through
a membrane module, so $\theta_{j}\in[0,1]$. However, when $\theta_{j}=1$,
$\ln(1-\theta_{j})$ (see (W18)) is not well-defined. Therefore, we choose
$\theta_{j}^{\textnormal{lo}}=0$ and
$\theta_{j}^{\textnormal{up}}=1-\epsilon_{\theta}$. For all our computations
in this article, we choose $\epsilon_{\theta}=10^{-3}$.

The choice of upper and lower bounds on mole fraction variables is listed in
(25). We choose the composition of the retentate product
($X^{\textnormal{out}}$) and the permeate product ($Y^{\textnormal{per}}$)
streams as the lower and upper bound on $x_{j}^{\textnormal{in}}$ (see (25a)),
respectively. This is justified, because each stage separates the
corresponding feed and a further separation would not be needed if the
composition of the feed is either above $Y^{\textnormal{per}}$ or below
$X^{\textnormal{out}}$. Next, we determine the bounds on
$y_{j}^{\textnormal{in}}$ using (8) and the bounds on
$x_{j}^{\textnormal{in}}$ and $u$. Since $y_{j}^{\textnormal{in}}$ increases
monotonically with $x_{j}^{\textnormal{in}}$ (see property P3 in §3.3) and $u$
(see §3.2), we obtain $(y_{j}^{\textnormal{in}})^{\textnormal{lo}}$ (resp.
$(y_{j}^{\textnormal{in}})^{\textnormal{up}}$) by substituting
$x^{\textnormal{in}}=(x_{j}^{\textnormal{in}})^{\textnormal{lo}}$ and
$u=u^{\textnormal{lo}}$ (resp.
$x^{\textnormal{in}}=(x_{j}^{\textnormal{in}})^{\textnormal{up}}$ and
$u=u^{\textnormal{up}}$) in (8) and solving for $y^{\textnormal{in}}$. This
procedure is symbolically represented as
$y|_{x=(x^{\textnormal{in}})^{\textnormal{lo}},\;u=u^{\textnormal{lo}}}$ in
(25b). Next, from property P5 in §4.2, the mole fraction of component $A$ in
the retentate decreases from stage 1 through $N$ i.e.,
$x_{N}^{\textnormal{out}}\leq\dots\leq x_{1}^{\textnormal{out}}$. Further,
since the retentate product is formed by mixing the retentate streams from
stages $N-1$ and $N$, its composition lies in the interval
$x^{\textnormal{out}}_{N}\leq X^{\textnormal{out}}\leq
x^{\textnormal{out}}_{N-1}$. Therefore, we impose $X^{\textnormal{out}}$ as
the lower bound on $x^{\textnormal{out}}_{j}$ for every
$j\in\\{1,\dots,N-1\\}$. While zero is a valid lower bound for
$x^{\textnormal{out}}_{N}$, specifying it leads to the following issue. When
$x^{\textnormal{out}}_{N}=0$, $y^{\textnormal{out}}_{N}=0$ from the flux
equation in (W18), and the term $\ln(z^{\textnormal{out}}_{N})$ in (W18) is
not well-defined. Therefore, we set
$(x^{\textnormal{out}}_{N})^{\textnormal{lo}}=\epsilon_{x}$ where
$\epsilon_{x}>0$ (see (25c)). For all our computations, we choose
$\epsilon_{x}=10^{-3}$. Next, as before, we obtain
$(y_{j}^{\textnormal{out}})^{\textnormal{lo}}$ using (9) and the bounds on
$x_{j}^{\textnormal{out}}$ and $u$ (see (25d)). Next, from property P4 in
§3.3, we have $y_{j}^{\textnormal{out}}\leq y_{j}^{\textnormal{per}}$. Thus,
the lower bound on $y_{j}^{\textnormal{out}}$ is a valid lower bound on
$y_{j}^{\textnormal{per}}$, so we choose
$(y_{j}^{\textnormal{per}})^{\textnormal{lo}}=(y_{j}^{\textnormal{out}})^{\textnormal{lo}}$
(see (25d)). Along the same line, the upper bound on
$y_{j}^{\textnormal{per}}$ is a valid upper bound on
$y_{j}^{\textnormal{out}}$, so we choose
$(y_{j}^{\textnormal{out}})^{\textnormal{up}}=(y_{j}^{\textnormal{per}})^{\textnormal{up}}$
(see (25e)). Further, $y_{j}^{\textnormal{per}}\leq y_{j}^{\textnormal{in}}$
from property P4 in §3.3, so the upper bound on $y_{j}^{\textnormal{in}}$ is a
valid upper bound on $y_{j}^{\textnormal{per}}$. However, for every
$j\in\\{2,\dots,N\\}$, a tighter upper bound can be inferred on
$y_{j}^{\textnormal{per}}$ from the following argument. Since
$y_{1}^{\textnormal{per}}\geq y_{2}^{\textnormal{per}}$ (see property P5 in
§4.2), $y_{2}^{\textnormal{per}}$ can be at most $Y^{\textnormal{per}}$ in
order to maintain the composition of the permeate product stream at
$Y^{\textnormal{per}}$. Thus, $Y^{\textnormal{per}}$ is a valid upper bound on
$y_{2}^{\textnormal{per}}$. Further, $Y^{\textnormal{per}}$ is also a valid
upper bound on $y_{3}^{\textnormal{per}},\dots,y_{N}^{\textnormal{per}}$
because $y_{2}^{\textnormal{per}}\geq\dots\geq y_{N}^{\textnormal{per}}$. This
leads to (25e). Finally, we obtain the upper bound on
$x_{j}^{\textnormal{out}}$ from (9) and the bounds on
$y_{j}^{\textnormal{out}}$ and $u=u^{\textnormal{lo}}$ as shown in (25f).

$\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle(x_{j}^{\textnormal{in}})^{\textnormal{lo}}=X^{\textnormal{out}},\qquad(x_{j}^{\textnormal{in}})^{\textnormal{up}}=Y^{\textnormal{per}}$
(25a) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle(y_{j}^{\textnormal{in}})^{\textnormal{lo}}=y|_{x=(x_{j}^{\textnormal{in}})^{\textnormal{lo}},\;u=u^{\textnormal{lo}}},\qquad(y_{j}^{\textnormal{in}})^{\textnormal{up}}=y|_{x=(x_{j}^{\textnormal{in}})^{\textnormal{up}},\;u=u^{\textnormal{up}}}$
(25b) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle(x_{j}^{\textnormal{out}})^{\textnormal{lo}}=\begin{cases}X^{\textnormal{out}},&\text{if
}j<N\\\ \epsilon_{x},&\text{otherwise}\end{cases}$ (25c)
$\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle(y_{j}^{\textnormal{per}})^{\textnormal{lo}}=(y_{j}^{\textnormal{out}})^{\textnormal{lo}}=y|_{(x_{j}^{\textnormal{out}})^{\textnormal{lo}},\;u=u^{\textnormal{lo}}}$
(25d) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle(y_{j}^{\textnormal{out}})^{\textnormal{up}}=(y_{j}^{\textnormal{per}})^{\textnormal{up}}=\begin{cases}(y_{j}^{\textnormal{in}})^{\textnormal{up}},&\text{if
}j=1\\\ Y^{\textnormal{per}},&\text{otherwise}\end{cases}$ (25e)
$\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle(x_{j}^{\textnormal{out}})^{\textnormal{up}}=x|_{y=(y_{j}^{\textnormal{out}})^{\textnormal{up}},\;u=u^{\textnormal{lo}}}$
(25f)

Lastly, we obtain the bounds on $z^{\textnormal{in}}_{j}$ and
$z_{j}^{\textnormal{out}}$ by analyzing the behavior of (8) and (9),
respectively. It can be verified that in the interval $[0,1]$, the RHS of both
the equations is concave, evaluates to zero at $y=0$ and $y=1$, and goes
through a maxima at $y=\frac{S-\sqrt{S}}{S-1}$. This leads to the choice of
bounds in (26). As before,
$z|_{y=(y_{j}^{\textnormal{in}})^{\textnormal{lo}},u=u^{\textnormal{lo}}}$
represents that the value of $z$ is obtained by substituting
$y=(y_{j}^{\textnormal{in}})^{\textnormal{lo}}$ and $u=u^{\textnormal{lo}}$ in
(8).

$\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle(z_{j}^{\textnormal{in}})^{\textnormal{lo}}=\min\left\\{z|_{y=(y_{j}^{\textnormal{in}})^{\textnormal{lo}},u=u^{\textnormal{lo}}},\;z|_{y=(y_{j}^{\textnormal{in}})^{\textnormal{up}},u=u^{\textnormal{lo}}}\right\\}$
(26a) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle(z_{j}^{\textnormal{out}})^{\textnormal{lo}}=\min\left\\{z|_{y=(y_{j}^{\textnormal{out}})^{\textnormal{lo}},u=u^{\textnormal{lo}}},\;z|_{y=(y_{j}^{\textnormal{out}})^{\textnormal{up}},u=u^{\textnormal{lo}}}\right\\}$
(26b) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle(z_{j}^{\textnormal{in}})^{\textnormal{up}}=\begin{cases}z|_{y=(y_{j}^{\textnormal{in}})^{\textnormal{up}},u=u^{\textnormal{up}}},&\text{if
}(y_{j}^{\textnormal{in}})^{\textnormal{up}}\leq\frac{S-\sqrt{S}}{S-1}\\\
k^{\textnormal{up}}(\sqrt{S}-1)^{2},&\text{if
}(y_{j}^{\textnormal{in}})^{\textnormal{lo}}\leq\frac{S-\sqrt{S}}{S-1}\leq(y_{j}^{\textnormal{in}})^{\textnormal{up}}\\\
z|_{y=(y_{j}^{\textnormal{in}})^{\textnormal{lo}},u=u^{\textnormal{up}}},&\text{if
}\frac{S-\sqrt{S}}{S-1}\leq(y_{j}^{\textnormal{in}})^{\textnormal{lo}}\end{cases}$
(26c) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle(z_{j}^{\textnormal{out}})^{\textnormal{up}}=\begin{cases}z|_{y=(y_{j}^{\textnormal{out}})^{\textnormal{up}},u=u^{\textnormal{up}}},&\text{if
}(y_{j}^{\textnormal{out}})^{\textnormal{up}}\leq\frac{S-\sqrt{S}}{S-1}\\\
k^{\textnormal{up}}(\sqrt{S}-1)^{2},&\text{if
}(y_{j}^{\textnormal{out}})^{\textnormal{lo}}\leq\frac{S-\sqrt{S}}{S-1}\leq(y_{j}^{\textnormal{out}})^{\textnormal{up}}\\\
z|_{y=(y_{j}^{\textnormal{out}})^{\textnormal{lo}},u=u^{\textnormal{up}}},&\text{if
}\frac{S-\sqrt{S}}{S-1}\leq(y_{j}^{\textnormal{out}})^{\textnormal{lo}}\end{cases}$
(26d)

###  4.5 Additional cuts

Here, we describe the additional constraints that are derived using the
properties P1 through P4 in §3.3 and P5 in §4.2. While some of these
constraints are redundant to the MINLP (W), they are not implied in the
relaxed problem where the nonlinear constraints appear in a relaxed form.
Providing these constraints explicitly helps global solvers in expediting the
convergence characteristics of the branch-and-bound algorithm. To the MINLP
(W), we append

$\displaystyle\forall\;j\in\mathcal{J},\quad$ $\displaystyle
x^{\textnormal{out}}_{j}\leq x^{\textnormal{in}}_{j}\leq
y_{j}^{\textnormal{per}},$ (W25) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle y_{j}^{\textnormal{out}}\leq y^{\textnormal{per}}_{j}\leq
y^{\textnormal{in}}_{j},$ (W26)
$\displaystyle\forall\;j\in\mathcal{J}\setminus\\{1\\},\quad$ $\displaystyle
x^{\textnormal{in}}_{j}\leq x^{\textnormal{in}}_{j-1},\quad
x^{\textnormal{out}}_{j}\leq x^{\textnormal{out}}_{j-1},\quad
y^{\textnormal{per}}_{j}\leq y^{\textnormal{per}}_{j-1}$ (W27)
$\displaystyle\forall\;j\in\mathcal{J}\setminus\\{1\\},\quad$ $\displaystyle
y^{\textnormal{in}}_{j}\leq y^{\textnormal{in}}_{j-1},\quad
y^{\textnormal{out}}_{j}\leq y^{\textnormal{out}}_{j-1}$ (W28) $\displaystyle
Y^{\textnormal{per}}\leq y^{\textnormal{per}}_{1}$ (W29) $\displaystyle
x^{\textnormal{out}}_{N}\leq X^{\textnormal{out}}$ (W30)

We obtain (i) (W25) using the properties P1 and P2, (ii) (W26) using the
property P4, (iii) (W27) using the property P5, and (iv) (W28) using (W27) and
the property P3. (W29) and (W30) are derived using the following arguments.
Observe that the permeate product stream is a mixture of the permeate streams
of stages 1 and 2. Since the upper bound on $y_{2}^{\textnormal{per}}$ is
$Y^{\textnormal{per}}$ (see (25e)), $y_{1}^{\textnormal{per}}$ must be at
least $Y^{\textnormal{per}}$ in order to maintain the composition of the
permeate product stream at $Y^{\textnormal{per}}$. Similarly, the retentate
product stream is a mixture of the retentate streams of stages $N-1$ and $N$.
Since the lower bound on $x_{N-1}^{\textnormal{out}}$ is
$X^{\textnormal{out}}$, $x_{N}^{\textnormal{out}}$ can be at most
$X^{\textnormal{out}}$ in order to maintain the composition of the retentate
product stream at $X^{\textnormal{out}}$. This concludes the formulation of
the MINLP.

###  4.6 Need for a Global Optimization Approach

Here, we demonstrate the need for a global optimization approach by showing
that local solvers can get trapped in suboptimal solutions even when discrete
variables are fixed to a specific cascade. As an example, consider the cascade
shown in Figure 7 to separate a liquid mixture of p-xylene and a
pseudocomponent (mixture of m-xylene and o-xylene). We choose the membrane
perm-selectivity to p-xylene to be 50. Let the molar flow rate and the
composition (mole fraction of p-xylene) of (i) the feed mixture be 250 mol/s
and 0.65, respectively (ii) the permeate product stream be 147 mol/s and
0.995. The molar flowrate and the composition of the retentate product stream
can be obtained from an overall mass balance. The values of the remaining
parameters are listed in the caption of Figure 7. Our objective is to identify
the optimal operating condition of the cascade that minimizes the overall
energy consumption. The MINLP (W) can be used for optimization after fixing
the binary variables to $\omega_{2}^{F}=1$,
$\omega^{\textnormal{recy},1}_{j}=1$ for $j=2,3,4$,
$\omega^{\textnormal{recy},2}_{j}=0$ for $j=3,4$,
$\omega^{\textnormal{per,bypass}}_{2}=0$, and
$\omega_{N-1}^{\textnormal{out,bypass}}=0$. When discrete variables are fixed,
(W) becomes a nonlinear program, so it can be solved using local solvers such
as CONOPT, SNOPT, IPOPT, etc. Here, we use CONOPT. To the best of our
knowledge, a systematic method is not available for identifying good initial
points. Therefore, we do not provide any initial point to the solver. With its
default initialization strategy, CONOPT yields an operating condition that
consumes 7,096 kW of power. Whereas, global solver BARON yields an operating
condition that consumes 1,780 kW of power ($\sim$75% reduction in power
consumption). This example clearly demonstrates the need to obtain the global
optimality certificate when optimizing a membrane cascade. Otherwise, we
cannot to be certain whether the obtained solution is the most energy
efficient or not.

Figure 7: Cascade considered in §4.6. The values of the missing parameters are
as follows: $V_{A}=1.233\times 10^{-4}$ m3/mol, $V_{B}=1.215\times 10^{-4}$
m3/mol, $\eta^{pump}=0.75$, $\eta^{TC}=0.80$, and $[(\Delta
P^{trans})^{\textnormal{lo}},(\Delta P^{trans})^{\textnormal{up}}]=$[30 bar,
107 bar]

###  4.7 Computational experiments

Here, through numerical experiments, we demonstrate that the proposed MINLP
(W) is able to identify the optimal membrane cascade within a relative
tolerance of 5%. In addition, we show the effectiveness of the cuts derived in
§4.5 in expediting the convergence characteristics of BARON by solving (W)
with and without (W25)–(W30). For our numerical experiments, we considered a
test set of 13 cases described in Table 3. The values of the remaining
parameters are reported in the caption. Note that instead of the flowrate and
the composition of permeate and retentate product streams, we have reported
the composition ($Y^{\textnormal{per}}$) and the recovery ($\gamma_{A}$) of
component $A$ in the permeate product stream. Recovery of a component is
defined as the ratio of its molar flowrate in the permeate product stream to
that in the feed. The values of the parameters needed in (W12)-(W15) can be
obtained from

$\displaystyle\gamma_{A}=\frac{F^{\textnormal{per}}\cdot
Y^{\textnormal{per}}}{F\cdot X^{F}},$ (27a) $\displaystyle
F=F^{\textnormal{per}}+F^{\textnormal{out}},$ (27b) $\displaystyle F\cdot
X^{F}=F^{\textnormal{per}}\cdot Y^{\textnormal{per}}+F^{\textnormal{out}}\cdot
X^{\textnormal{out}}.$ (27c)

(27a) is the definition of the recovery of component $A$, (27b) is the overall
mass balance across the superstructure, and (27c) is the mass balance of
component $A$ across the superstructure. For all gaseous (resp. liquid)
mixtures, we choose the admissible range of operating pressure ratio (resp.
trans-membrane pressure difference) to be $[1.1,9]$ (resp. $[30\text{
bar},107\text{ bar}]$) i.e., $r^{\textnormal{lo}}=1.1$ and
$r^{\textnormal{up}}=9$ (resp. $(\Delta P^{trans})^{\textnormal{lo}}=30$ bar
and $(\Delta P^{trans})^{\textnormal{up}}=107$ bar).

We use BARON 18.5.8 on GAMS 25.1 to solve the MINLP (W). All BARON options
except pDo were left at their default values. pDo was set to $-1$. We set the
relative tolerance for convergence ($\epsilon_{r}$) to 5% and the time limit
to 40 h as the termination criteria. All computations are performed on a Dell
Optiplex 5040 with 16 GB RAM, which has Intel Core i7-6700 3.4 GHz processor
and is running 64-bit Windows 10. The computational results are summarized in
Table 3.

The eighth (resp. ninth) column in Table 3 lists the computational performance
when the MINLP (W) is solved without (resp. with) (W25)–(W30). Clearly, when
the additional cuts are not included, none of the cases converge even after 40
hours. The remaining duality gap (defined as (Best known upper bound $-$ Best
known lower bound)/Best known upper bound) at the end of 40 hours is as high
as 98% in some cases. On the other hand, with the inclusion of (W25)–(W30), we
could solve all 13 cases to 5%-optimality within 40 hours. It is interesting
to note that without (W25)–(W30), the obtained solutions without optimality
certificate, were also at optimality as gleaned by comparing these solutions
with those obtained when appending (W25)–(W30). However, in order to have
confidence of weather the obtained solutions are globally optimal, it is
essential to obtain the optimality certificate within the desired optimality
gap. It avoids cases where one may have suboptimal solutions at hand.
Therefore, we recommend to always use constraints (W25)–(W30) when solving the
postulated cascade optimization problem.

Table 3: Test set for computational experiments. All mole fractions correspond to the most permeable component. The missing parameters are set to $N=4$, $F=250$ mol/s, $\eta^{comp}=0.75$, $\eta^{pump}=0.75$, $\eta^{TC}=0.8$, $V_{A}=1.233\times 10^{-4}$ m3/mol, $V_{B}=1.215\times 10^{-4}$ m3/mol, and $T=303.15$ K. Molar volumes are determined using Aspen Plus v8.6, and the reported $V_{B}$ is the average of molar volumes of o-xylene and m-xylene. All cases were solved with BARON 18.5.8 in GAMS 25.1 with a time limit of 40 hours. Except epsr and pDo, all other BARON options were left at their default values. We set epsr = 0.05 and pDo = -1. TL indicates that the time limit was reached, and the number in parenthesis is the remaining duality gap at the end of 40 hours. I/T correspond to number of BARON iterations/computational time (in hours) to reach 5%-optimality gap.  Case |  Mixture |  Feed |  Feed |  Permeate |  Molar |  Membrane |  Without |  With
---|---|---|---|---|---|---|---|---
|  |  Pressure |  Composition |  Purity |  Recovery |  perm- |  (W25)–(W30) |  (W25)–(W30)
|  |  (bar) |  (% mol) |  (% mol) |  (%) |  selectivity |  I/T |  I/T
Gaseous Mixtures
1 |  CO2/CH4 |  20 |  60 |  95 |  98.7 |  25 |  TL (75%) |  67,523 / 32.8
2 |  |  20 |  10 |  95 |  81.7 |  25 |  TL (89%) |  14,297 / 3.2
3 |  |  20 |  10 |  69 |  82.4 |  12 |  TL (92%) |  16,173 / 3.5
4 |  |  55 |  30 |  86.7 |  93.2 |  21 |  TL (74%) |  12,537 / 3.0
5 |  H2/CO2 |  28 |  21.5 |  99 |  98.2 |  38 |  TL (94%) |  26,627 / 4.8
6 |  propylene / propane |  9 |  80 |  90 |  90 |  5 |  TL (70%) |  51,845 / 16.6
7 |  |  9 |  70 |  99.6 |  97.8 |  35 |  TL (89%) |  66,903 / 20.8
8 |  |  9 |  70 |  92 |  97.8 |  35 |  TL (62%) |  33,211 / 13.2
Liquid Mixtures
9 |  p-xylene / pseudocomponent |  1 |  65 |  99.5 |  90 |  50 |  TL (72%) |  6,631 / 2.2
10 |  |  1 |  65 |  |  99 |  |  TL (93%) |  6,963 / 1.8
11 |  |  1 |  90 |  |  90 |  |  TL (64%) |  25,848 / 14.5
12 |  |  1 |  90 |  |  99 |  |  TL (72%) |  9,091 / 2.7
13 |  |  1 |  23.6 |  |  97.5 |  |  TL (98%) |  4,331 / 1.0

##  5 Case studies

Here, we examine Case 8 (separation of propylene/propane mixture) and Case 12
(separation of p-xylene/(o-xylene+m-xylene) mixture) in Table 3 in more
detail. In each case, we determine the optimal membrane cascade requiring at
most four stages along with its optimal operating conditions using the MINLP
(W). Additionally, we derive constraints that can be appended to the MINLP to
limit the number of compressors/pumps in the optimal cascade. All optimization
problems in this section are solved to 5%-optimality gap.

###  5.1 Case 8: Propylene/propane separation

The separation of propylene/propane mixture is of great relevance at
industrial level, as propylene is an important precursor in the production of
a wide variety of chemicals such as polypropylene, propylene oxide, acrylic
acid, acrylonitrile, and isopropanol [52]. Currently, a sub-ambient
distillation at mild pressure is used to separate propylene/propane mixtures
[53, 54]. Given the very high demand of propylene, it is typically separated
with a high recovery that ranges between 97.8% to 99% [55, 56]. On the other
hand, the final purity of the propylene product varies according to its use in
the downstream processes. Propylene needed for the synthesis of polymers is
generally produced at $>$99.5% purity. Whereas, propylene needed for the
synthesis of chemicals, such as 2-propanol, is produced at moderate purity of
about 92% [52, 57].

For our case study, we consider the separation of a 70%/30% mixture of
propylene/propane to produce chemical grade propylene at 92% purity and 97.8%
recovery. The chosen feed composition and product recovery were taken from
[55], and they are in good agreement with the values encountered in industry
[58]. We take the molar flow rate of the feed to be 250 mol/s, and we assume
that the feed is available at 9 bar. Further, we assume that all the membranes
in the cascade have a perm-selectivity of 35 [55] (note that propylene is the
most permeable component).

####  5.1.1 Optimal four-stage cascade

The optimal four-stage membrane cascade for the separation of the
propylene/propane mixture is shown in Figure 8. In this cascade, the feed
mixture is located at the first stage. This is not surprising given that the
degree of enrichment of the most permeable component in the permeate product
stream is not too high. Thus, a single stage is sufficient to achieve the
desired purity specification. Nevertheless, it is not sufficient to recover
the desired amount of propylene while maintaining the purity. As shown in
Figure 8, the optimal cascade uses three additional stages for stripping
propylene from the retentate stream leaving the first stage. At optimal
operation, the purity of the permeate from the first stage is higher than the
desired purity (see Figure 8). Thus, the permeate from the second stage, which
does not satisfy the purity requirement (see Figure 8), can be mixed with that
from the first stage. This enables the production of the permeate product
stream at the desired purity while avoiding the recycling of the permeate from
the second stage.

Figure 8: Optimal four-stage cascade for the separation of the
propylene/propane mixture. All mole fractions shown in the box-legends
correspond to propylene. Any discrepancy in the mass balances is due to the
rounding of the flow and the composition values.

####  5.1.2 Optimal cascade with only one intermediate compressor

The optimal four-stage cascade has two intermediate compressors, and one
compressor to pressurize the permeate product stream. Owing to the high
capital cost and maintenance issues associated with compressors, process
designers often seek cascades with few compressors. In the following, we
formulate the constraints which can be appended to the MINLP (W) to obtain
cascades with a specified maximum number of compressors.

Excluding the compressor needed to compress the permeate product stream, the
maximum number of compressors, $M_{comp}$, in a cascade can be controlled
using the binary variables $\omega_{j}^{\textnormal{recy},1}$ and
$\omega_{j}^{\textnormal{recy},2}$ through

$\displaystyle\sum_{j=2}^{N}\omega_{j}^{\textnormal{recy},1}+\sum_{j=3}^{N}\omega_{j}^{\textnormal{recy},2}-\sum_{j=3}^{N}\omega_{j-1}^{\textnormal{recy},1}\omega_{j}^{\textnormal{recy},2}\leq
M_{comp},$ (28)

The LHS computes the number of compressors in the cascade. The first (resp.
second) term in the above equation accounts the compressor along the arc
connecting the splitter $\mathcal{S}_{j}$ and mixer $\mathcal{M}_{j-1}$ (resp.
$\mathcal{M}_{j-2}$) (see Figure 5(a)). In cases where the material flow from
both $\mathcal{S}_{j-1}$ and $\mathcal{S}_{j}$ is sent to mixer
$\mathcal{M}_{j-2}$, both the streams can be mixed before compressing. This
reduces the number of compressors needed for the cascade. Thus, the third term
in the above equation reduces the number of compressors by one, when both
$\omega_{j-1}^{\textnormal{recy},1}=1$ and
$\omega_{j}^{\textnormal{recy},2}=1$. We linearize (28) by introducing an
auxiliary variable $h_{j-1}$ for the product, and relax
$h_{j-1}=\omega_{j-1}^{\textnormal{recy},1}\omega_{j}^{\textnormal{recy},2}$
by replacing the bilinear term with its convex and concave envelopes over
$[0,1]^{2}$. This results in the following constraints

$\displaystyle\sum_{j=1}^{N}\omega_{j}^{\textnormal{recy},1}+\sum_{j=3}^{N}\omega_{j}^{\textnormal{recy},2}-\sum_{j=3}^{N}h_{j-1}\leq
M_{comp},$ (W31)
$\displaystyle\forall\;j\in\left\\{3,\dots,N\right\\},\quad\left\\{\begin{aligned}
&h_{j-1}\leq\omega_{j-1}^{\textnormal{recy},1},\quad
h_{j-1}\leq\omega_{j}^{\textnormal{recy},2}\\\ &h_{j-1}\geq 0,\quad
h_{j-1}\geq\omega_{j-1}^{\textnormal{recy},1}+\omega_{j}^{\textnormal{recy},2}-1,\end{aligned}\right.$
(W32)

Since $\omega_{j-1}^{\textnormal{recy},1}$ and
$\omega_{j}^{\textnormal{recy},2}$ are binary, it can be easily verified that
the feasible solutions to (W31) and (W32) are the same as that of (28).

Suppose that a stage $j$ does not carryout any separation i.e.,its stage cut
is zero ($\theta_{j}=0$). Then, the molar flow rate of the permeate stream
from stage $j$ is zero, so the compressors along the arcs connecting the
splitter $\mathcal{S}_{j}$ and the mixers $\mathcal{M}_{j-1}$ and
$\mathcal{M}_{j-2}$ are not needed. However, (W21) does not permit the
optimizer to set both $\omega_{j}^{\textnormal{recy},1}=0$ and
$\omega_{j}^{\textnormal{recy},2}=0$. Instead it requires that either
$\omega_{j}^{\textnormal{recy},1}=1$ or $\omega_{j}^{\textnormal{recy},2}=1$,
which contributes towards the number of compressors in the cascade. Therefore,
we modify (W21) as shown in (W33). Here, $\mu_{j}$ is a binary variable which
takes the value one when the stage cut of a stage is nonzero (see (W34)). This
way, the optimizer can suppress the compressors when the corresponding stage
cut becomes zero by setting $\mu_{j}=0$.

$\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle\left.\begin{aligned}
&\omega^{\textnormal{recy},1}_{j}\delta_{j\geq
2}+\omega^{\textnormal{recy},2}_{j}\delta_{j\geq
3}+\omega^{\textnormal{per,bypass}}_{j}\delta_{j\leq 2}=\mu_{j},\\\
&\omega_{j}^{\textnormal{recy},1}\in\\{0,1\\},\quad\omega_{j}^{\textnormal{recy},2}\in\\{0,1\\},\quad\omega_{j}^{\textnormal{per,bypass}}\in\\{0,1\\},\quad\mu_{j}\in\\{0,1\\}\end{aligned}\right\\}$
(W33) $\displaystyle\forall\;j\in\mathcal{J},\quad$
$\displaystyle\theta_{j}\leq\mu_{j},$ (W34)

We now append (W31)-(W34) to (W), eliminate (W21), and solve the resulting
MINLP by substituting $M_{comp}=1$ to obtain the optimal cascade requiring
only one intermediate compressor. The optimal solution, shown in Figure 9,
requires only three stages. The optimizer bypasses material flow through a
stage by setting its stage cut to zero. This cascade requires only 1.3% more
power than the one in Figure 8. Since the capital cost of a compressor can be
significant, the small increase in energy cost will be dominated by the
reduction in the overall capital cost making this arrangement more economical
compared to the one in Figure 8. A quantitative cost benefit analysis of the
cascades in Figures 8 and 9 is beyond the scope of this paper.

Figure 9: Optimal cascade with only one intermediate compressor for the
separation of the propylene/propane mixture. All mole fractions shown in the
box-legends correspond to propylene. Any discrepancy in the mass balances is
due to the rounding of the flow and the composition values.

###  5.2 Case 12: Separation of a mixture of xylenes

Xylenes are important compounds used in the manufacturing of a variety of
products ranging from plastics, paint solvents, resins, inks, and plasticizers
[59]. Among all the isomers, p-xylene has the highest commercial demand. One
of the major applications of p-xylene is the synthesis of terephthalic acid,
which is an important building block for manufacturing polymers such as
polyethylene terephthalate (PET) and polybutylene terephthalate (PBT) [60].
p-Xylene is primarily produced by catalytic reforming of naphtha. In addition,
it is also produced via the toluene disproportionation process. However,
compared to catalytic reforming, the latter process is more attractive from
the following perspective. The C8 cut obtained from the reforming process
contains only 20 mol% of p-xylene. The rest of the mixture consists of
o-xylene, m-xylene and ethyl benzene (EB). On the other hand, the C8 cut
obtained from the toluene disproportionation process contains around 90 mol%
of p-xylene[61, 62], and EB is practically absent in the cut [59]. In either
case, a separation is needed to concentrate p-xylene since it is needed at
$>99\%$ purity for the downstream processes.

Currently, simulated moving bed and crystallization are predominantly used to
separate p-xylene from the C8 cut. The inherent energy intensive nature of the
separation is due to similarities in physical and structural properties of the
constituent components, and it motivates the development of alternative
energy-efficient technologies. Recently, reverse osmosis is proposed as a
plausible alternative [63], so we investigate the separation of p-xylene from
the C8 cut using a membrane cascade. We assume that the feed mixture comes
from the toluene disproportionation process, so the composition of ethyl
benzene in the mixture is negligible. Further, we treat the C8 cut as a binary
mixture by lumping both o-xylene and m-xylene as a single pseudo component.
This is justified, because the membranes available in the literature are more
selective to p-xylene compared to both o-xylene and m-xylene [63, 64].

For our case study, we consider a 250 mol/s mixture of xylenes containing 90%
p-xylene at 30 °C and 1 bar. Since p-xylene is a valuable chemical and it is
needed at a high purity for downstream processes, we require that the permeate
product stream recovers 99% of p-xylene at 99.5% purity [62]. The perm-
selectivity of p-xylene w.r.t o-xylene+m-xylene pseudo component is taken to
be 50. This value is close to the experimentally observed selectivity of
p-xylene w.r.t o-xylene through a carbon molecular sieve membrane (see Koh et
al. [63]). Finally, we let the admissible operating range of the trans-
membrane pressure difference be 30 bar – 107 bar.

####  5.2.1 Optimal four-stage cascade

Figure 10 shows the optimal four-stage cascade along with its optimal
operating condition. In this cascade, the feed mixture is fed to the second
stage. With the exception of the first stage, the permeate stream from each
stage is pumped and recycled to the stage on its left. For this reason, the
cascade requires three pumps to pressurize each of the recycle streams. This
is in addition to the pump and the turbocharger needed to pressurize the feed
mixture (see Figure 10).

Figure 10: Optimal four-stage cascade for the separation of xylene isomers.
All mole fractions shown in the box-legends correspond to p-xylene. Any
discrepancy in the mass balances is due to the rounding of the flow and the
composition values.

####  5.2.2 On the admissible operating range of $\Delta P^{trans}$

In Figure 10, at optimal operation, we observe that the trans-membrane
pressure difference ($\Delta P^{trans}$) is at its upper bound. This suggests
that the overall energy consumption could be reduced by increasing the upper
bound on $\Delta P^{trans}$. To verify this hypothesis, we increase
$u^{\textnormal{up}}$ and solve the MINLP. The results are shown in Table 4.
As expected, the total power requirement reduces when the cascade is allowed
to operate at a higher $\Delta P^{trans}$. Unfortunately, the membranes that
are currently available suffer from compaction and lose selectivity towards
p-xylene dramatically approximately beyond 107 bar111In the experimental
permeation data provided by Prof. Ryan Lively, we observed that membrane perm-
selectivity drops when the trans-membrane pressure difference exceeds 107
bar.[63]. Therefore, novel membrane materials with improved mechanical
properties are needed to further reduce the power requirement.

$u^{\textnormal{up}}$ (bar) |  Optimum $\Delta P^{trans}$ (bar) |  Total power (kW)
---|---|---
107 |  107 |  1,416
125 |  125 |  1,281
150 |  150 |  1,205
Table 4: Optimal power consumption for the cascade in Figure 10 for different
ranges of admissible values of $\Delta p^{trans}$.

####  5.2.3 Determining target membrane selectivity

An alternative, and well-known, approach to reduce the energy requirement of a
cascade is to increase the perm-selectivity of the membrane. Currently, the
membrane literature relies heavily on the Robeson plot to assess the
performance of a membrane [65]. Unfortunately, this analysis does not answer
how the energy requirement of a membrane cascade compares with an alternate
technology. Especially, if membranes are being developed in order to replace
an existing technology, it is useful to determine the minimum perm-selectivity
needed to compete with the existing alternate technology in terms of the
energy consumption. We now show how to determine the target selectivity for a
given application using our MINLP formulation.

We solve the MINLP (W) for different values of selectivity and obtain the
optimal power requirement for each value. The results are summarized in Table
5. For the separation of a mixture of xylenes, the state-of-the-art polymeric
membranes available in the literature have a selectivity of $\sim$50\. For
such membranes, the optimal power requirement is 1,416 kW. If the existing
non-membrane separation technology consumes more power, then replacing it with
a membrane cascade would be beneficial from the point of view of energy
consumption. On the other hand, if the existing separation technology consumes
less power, say by 10%, then replacing it with a membrane cascade would be
detrimental. In that case, a membrane perm-selectivity of at least 63 is
needed to compete with the existing technology. Similarly, if the existing
technology requires 20% (resp. 40%) less power, then a membrane perm-
selectivity of at least 81 (resp. 231) is needed in order for the membrane
cascade to operate with the same power as the existing technology. We note
that energy is only one component of the overall separation cost. While we
only discuss energy here, in practice, one will have to additionally account
for the capital costs for a proper comparison.

p-xylene perm-selectivity |  Total power (kW)
---|---
50 |  1,416
63 |  1,274
81 |  1,133
231 |  850
Table 5: Optimum power consumption for the cascade in Figure 10 for different
values of p-xylene perm-selectivity.

####  5.2.4 Optimal four-stage cascade with two intermediate pumps

Unlike compressors, pumps are relatively cheap and easy to operate. Thus, the
number of pumps used for liquid separations may not be a major concern.
Nevertheless, operating a cascade beyond a certain number of pumps may not be
economically attractive. In such cases, one can append (W31)-(W34) to the
MINLP and solve it to obtain attractive cascades with limited number of pumps.

For the xylene mixture, we determine the optimal cascade requiring at most two
intermediate pumps i.e., $M_{comp}=2$. Solving (W) after appending (W31)-(W34)
yields the solution shown in Figure 11. Similar to Figure 10, the feed mixture
is fed at the second stage in Figure 11. However, in the latter arrangement,
the permeate streams from both the second and third stages are mixed and
recycled to the first stage. At optimum operation, this cascade requires 13%
more energy than the one in Figure 10. An economic analysis of these
arrangements reveals whether the reduction in the capital cost due to
elimination of a pump overweights the increase in operating cost due to the
increase in power consumption or not.

Figure 11: Optimal four-stage cascade with only two intermediate pumps for the
separation of xylene isomers. All mole fractions shown in the box-legends
correspond to p-xylene. Any discrepancy in the mass balances is due to the
rounding of the flow and the composition values.

####  5.2.5 Optimal three-stage cascade

The main purpose of this subsection is to demonstrate that solving the MINLP
with (W31)-(W34) often leads to better solutions than simply solving (W) with
a fewer number of stages. For the xylene separation, a membrane cascade with
only two intermediate pumps can also be designed by solving (W) with $N=3$.
Figure 12 shows the optimal cascade along with its optimal operating condition
obtained by solving (W) with $N=3$. Interestingly, the optimal three-stage
cascade consumes significantly more energy than the cascade in Figure 11 even
though both the arrangements require the same number of pumps. This example
illustrates that (W31)-(W34) allow discovery of energy efficient cascades with
limited number of pumps/compressors which are otherwise not easy to obtain
directly by solving the MINLP without these constraints.

Figure 12: Optimal three-stage cascade for the separation of p-xylene. All
mole fractions shown in the box-legends refer to p-xylene. Any discrepancy in
the mass balances is due to the rounding of the flow and the composition
values.

##  6 Concluding remarks

In this article, we addressed the problem of identifying energy efficient
membrane cascades for the separation of a given binary liquid or gaseous
mixture through dense membranes. The main highlights of the article are the
following. First, to model the permeation process, we used a simple yet
reasonably accurate model that does not require the solution of a detailed set
of transport equations on either side of the membrane. This model was derived
by solving a DAE system analytically. An advantage of the analytical solution
is that it circumvents the need for discretization, where the latter technique
requires a trade-off between model accuracy and complexity. Second, we
proposed a novel MINLP that is applicable for both gaseous and liquid
mixtures. Third, we demonstrated the need for a global optimality certificate
on a numerical example, where we showed that local solution may consume
significantly more energy than the globally optimal cascade. Unfortunately,
state-of-the-art global solvers fail to solve the problem to a specified
optimality tolerance in a reasonable amount of time. To aid the convergence of
the solvers, we derived additional cuts that utilize physical intuition and
the mathematical structure of the equations. Our computational experiments
demonstrated that these cuts significantly improve the convergence of the
global solvers. Finally, we used our model to find attractive cascades for two
important industrial separations.

## Acknowledgements

Jose Adrian Chavez Velasco thanks COMEXUS-Fulbright, and CONACYT for providing
him with financial support. We thank Professor Ryan Lively for sharing data in
reference [63] with us.

##  Appendix A Appendix

###  A.1 Reformulation of the flux equation

We begin with the definition of the local mole fraction on the permeate side

$\displaystyle
y=\frac{n_{A}}{n_{A}+n_{B}}=\frac{S\left[x-y\exp(-C_{A}u)\right]}{S\left[x-y\exp(-C_{A}u)\right]+\left[(1-x)-(1-y)\exp(-C_{B}u)\right]}$
$\displaystyle\implies$ $\displaystyle y=\frac{Sx-
Sy\exp(-C_{A}u)}{x(S-1)-y\left[S\exp(-C_{A}u)-\exp(-C_{B}u)\right]+1-\exp(-C_{B}u)}$

Multiplying both sides with the denominator of the right hand side leads to

$\displaystyle yx(S-1)+y=Sx-\left[S\exp(-C_{A}u)-\exp(-C_{B}u)\right]y(1-y)$
$\displaystyle\implies$
$\displaystyle\frac{y}{S-(S-1)y}=x-\frac{\left[S\exp(-C_{A}u)-\exp(-C_{B}u)\right]y(1-y)}{S-(S-1)y}$
$\displaystyle\implies$ $\displaystyle
y-x=y-\frac{y}{S-(S-1)y}-\frac{\left[S\exp(-C_{A}u)-\exp(-C_{B}u)\right]y(1-y)}{S-(S-1)y}$
$\displaystyle\implies$ $\displaystyle
y-x=\left[S(1-e^{-C_{A}u})-(1-e^{-C_{B}u})\right]\frac{y(1-y)}{S-(S-1)y}$

we use partial fraction decomposition to rearrange the right hand side as
follows

$\displaystyle
y-x=\left[\frac{S(1-e^{-C_{A}u})-(1-e^{-C_{B}u})}{(S-1)^{2}}\right]\left(1+(S-1)y-\frac{S}{S-(S-1)y}\right).$

Finally, we obtain (5) by defining the term inside square brackets as $k$.

###  A.2 Derivation of the Analytical Solution

Consider the DAE system

$\displaystyle\frac{dx}{df}=\frac{y-x}{f},\quad
x(f^{\textnormal{in}})=x^{\textnormal{in}},\quad
f\in[f^{\textnormal{in}},f^{\textnormal{out}}]$ $\displaystyle
y-x=k\left(1+(S-1)y-\frac{S}{S-(S-1)y}\right)$

This implies

$\displaystyle\int_{x^{\textnormal{in}}}^{x^{\textnormal{out}}}\frac{dx}{y-x}=\int_{x^{\textnormal{in}}}^{x^{\textnormal{out}}}\frac{d(y-(y-x))}{y-x}=\int_{f^{\textnormal{in}}}^{f^{\textnormal{out}}}\frac{df}{f},$
$\displaystyle\implies$
$\displaystyle\int_{y^{\textnormal{in}}}^{y^{\textnormal{out}}}\frac{dy}{y-x}-\int_{y^{\textnormal{in}}-x^{\textnormal{in}}}^{y^{\textnormal{out}}-x^{\textnormal{out}}}\frac{d(y-x)}{y-x}=\int_{f^{\textnormal{in}}}^{f^{\textnormal{out}}}\frac{df}{f}$
$\displaystyle\implies$
$\displaystyle\frac{1}{k(S-1)^{2}}\int_{y^{\textnormal{in}}}^{y^{\textnormal{out}}}\frac{S-(S-1)y}{y(1-y)}dy-\int_{y^{\textnormal{in}}-x^{\textnormal{in}}}^{y^{\textnormal{out}}-x^{\textnormal{out}}}\frac{d(y-x)}{y-x}=\int_{f^{\textnormal{in}}}^{f^{\textnormal{out}}}\frac{df}{f}$
$\displaystyle\implies$
$\displaystyle\frac{1}{k(S-1)^{2}}\int_{y^{\textnormal{in}}}^{y^{\textnormal{out}}}\left[\frac{S}{y}+\frac{1}{1-y}\right]dy-\ln\frac{y^{\textnormal{out}}-x^{\textnormal{out}}}{y^{\textnormal{in}}-x^{\textnormal{in}}}=\ln\frac{f^{\textnormal{out}}}{f^{\textnormal{in}}}$
$\displaystyle\implies$
$\displaystyle\frac{1}{k(S-1)^{2}}\left[S\ln\frac{y^{\textnormal{out}}}{y^{\textnormal{in}}}-\ln\frac{1-y^{\textnormal{out}}}{1-y^{\textnormal{in}}}\right]-\ln\frac{y^{\textnormal{out}}-x^{\textnormal{out}}}{y^{\textnormal{in}}-x^{\textnormal{in}}}=\ln\frac{f^{\textnormal{out}}}{f^{\textnormal{in}}}$

###  A.3 Derivation of $dy/dx$ used in Property P3

We differentiate (5) with respect to $y$ to obtain

$\displaystyle\frac{dx}{dy}$
$\displaystyle=1-k(S-1)+\frac{k(S-1)S}{\left[S-(S-1)y\right]^{2}}$
$\displaystyle=\frac{1}{y}\left\\{y-k(S-1)y+\frac{ky(S-1)S}{\left[S-(S-1)y\right]^{2}}\right\\}$
$\displaystyle\overset{\eqref{genfluxrel_2}}{=}\frac{1}{y}\left\\{x+k-\frac{kS}{S-(S-1)y}+\frac{ky(S-1)S}{\left[S-(S-1)y\right]^{2}}\right\\}$
$\displaystyle=\frac{1}{y}\left\\{x+\frac{k(S-1)^{2}y^{2}}{\left[S(1-y)+y\right]^{2}}\right\\}.$

Therefore,

$\displaystyle\frac{dy}{dx}=\left(\frac{dx}{dy}\right)^{-1}=y\left[x+\frac{k(S-1)^{2}y^{2}}{\left[S(1-y)+y\right]^{2}}\right]^{-1}.$

## References

  * [1] M. Galizia, W. S. Chi, Z. P. Smith, T. C. Merkel, R. W. Baker, and B. D. Freeman, “50th anniversary perspective: Polymers and mixed matrix membranes for gas and vapor separation: A review and prospective opportunities,” Macromolecules, vol. 50, no. 20, pp. 7809–7843, 2017.
  * [2] L. S. White, “Development of large-scale applications in organic solvent nanofiltration and pervaporation for chemical and refining processes,” Journal of Membrane Science, vol. 286, no. 1, pp. 26 – 35, 2006.
  * [3] R. Spillman, “Chapter 13 economics of gas separation membrane processes,” vol. 2 of Membrane Science and Technology, pp. 589 – 667, Elsevier, 1995\.
  * [4] R. Agrawal, “A simplified method for the synthesis of gas separation membrane cascades with limited numbers of compressors,” Chemical Engineering Science, vol. 52, no. 6, pp. 1029 – 1044, 1997.
  * [5] R. Qi and M. Henson, “Optimization-based design of spiral-wound membrane systems for CO2/CH4 separations,” Separation and Purification Technology, vol. 13, no. 3, pp. 209 – 225, 1998\.
  * [6] J. A. Caballero and I. E. Grossmann, “Structural considerations and modeling in the synthesis of heat-integrated-thermally coupled distillation sequences,” Industrial & Engineering Chemistry Research, vol. 45, no. 25, pp. 8454–8474, 2006.
  * [7] Z. Jiang, G. Madenoor Ramapriya, M. Tawarmalani, and R. Agrawal, “Minimum energy of multicomponent distillation systems using minimum additional heat and mass integration sections,” AIChE Journal, vol. 64, no. 9, pp. 3410–3418, 2018.
  * [8] U. Nallasivam, V. H. Shah, A. A. Shenvi, J. Huff, M. Tawarmalani, and R. Agrawal, “Global optimization of multicomponent distillation configurations: 2. Enumeration based global minimization algorithm,” AIChE Journal, vol. 62, no. 6, pp. 2071–2086, 2016.
  * [9] R. T. Gooty, R. Agrawal, and M. Tawarmalani, “An MINLP formulation for the optimization of multicomponent distillation configurations,” Computers & Chemical Engineering, vol. 125, pp. 13 – 30, 2019.
  * [10] M. Tawarmalani and N. V. Sahinidis, “A polyhedral branch-and-cut approach to global optimization,” Mathematical programming, vol. 103, no. 2, pp. 225–249, 2005.
  * [11] M. R. Kılınç and N. V. Sahinidis, “Exploiting integrality in the global optimization of mixed-integer nonlinear programming problems with baron,” Optimization Methods and Software, vol. 33, no. 3, pp. 540–562, 2018.
  * [12] M. M. Qiu, S. T. Hwang, and Y. K. Kao, “Economic evaluation of gas membrane separator designs,” Industrial & Engineering Chemistry Research, vol. 28, no. 11, pp. 1670–1677, 1989.
  * [13] B. Bhide and S. Stern, “A new evaluation of membrane processes for the oxygen-enrichment of air. I. Identification of optimum operating conditions and process configuration,” Journal of Membrane Science, vol. 62, no. 1, pp. 13 – 35, 1991.
  * [14] B. Bhide and S. Stern, “Membrane processes for the removal of acid gases from natural gas. I. Process configurations and optimization of operating conditions,” Journal of Membrane Science, vol. 81, no. 3, pp. 209 – 237, 1993.
  * [15] J. Xu and R. Agrawal, “Gas separation membrane cascades I. One-compressor cascades with minimal exergy losses due to mixing,” Journal of Membrane Science, vol. 112, no. 2, pp. 115 – 128, 1996.
  * [16] R. Agrawal and J. Xu, “Gas separation membrane cascades II. Two-compressor cascades,” Journal of Membrane Science, vol. 112, no. 2, pp. 129 – 146, 1996.
  * [17] H. Lababidi, G. A. Al-Enezi, and H. M. Ettouney, “Optimization of module configuration in membrane gas separation,” Journal of Membrane Science, vol. 112, no. 2, pp. 185 – 197, 1996.
  * [18] J. Hao, P. Rice, and S. Stern, “Upgrading low-quality natural gas with H2S- and CO2-selective polymer membranes: Part II. Process design, economics, and sensitivity study of membrane stages with recycle streams,” Journal of Membrane Science, vol. 320, no. 1, pp. 108 – 122, 2008.
  * [19] F. Ahmad, K. Lau, A. Shariff, and G. Murshid, “Process simulation and optimal design of membrane separation system for CO2 capture from natural gas,” Computers & Chemical Engineering, vol. 36, pp. 119 – 128, 2012.
  * [20] R. Agrawal and J. Xu, “Gas-separation membrane cascades utilizing limited numbers of compressors,” AIChE Journal, vol. 42, no. 8, pp. 2141–2154, 1996.
  * [21] R. Pathare and R. Agrawal, “Design of membrane cascades for gas separation,” Journal of Membrane Science, vol. 364, no. 1, pp. 263 – 277, 2010.
  * [22] A. Aliaga-Vicente, J. A. Caballero, and M. J. Fernández-Torres, “Synthesis and optimization of membrane cascade for gas separation via mixed-integer nonlinear programming,” AIChE Journal, vol. 63, no. 6, pp. 1989–2006, 2017\.
  * [23] R. Qi and M. A. Henson, “Optimal design of spiral-wound membrane networks for gas separations,” Journal of Membrane Science, vol. 148, no. 1, pp. 71 – 89, 1998.
  * [24] R. Qi and M. A. Henson, “Membrane system design for multicomponent gas mixtures via mixed-integer nonlinear programming,” Computers & Chemical Engineering, vol. 24, no. 12, pp. 2719 – 2737, 2000.
  * [25] R. V. S. Uppaluri, P. Linke, and A. C. Kokossis, “Synthesis and optimization of gas permeation membrane networks,” Industrial & Engineering Chemistry Research, vol. 43, no. 15, pp. 4305–4322, 2004.
  * [26] R. V. Uppaluri, R. Smith, P. Linke, and A. C. Kokossis, “On the simultaneous optimization of pressure and layout for gas permeation membrane systems,” Journal of Membrane Science, vol. 280, no. 1, pp. 832 – 848, 2006.
  * [27] M. Scholz, M. Alders, T. Lohaus, and M. Wessling, “Structural optimization of membrane-based biogas upgrading processes,” Journal of Membrane Science, vol. 474, pp. 1 – 10, 2015.
  * [28] Álvaro A. Ramírez-Santos, M. Bozorg, B. Addis, V. Piccialli, C. Castel, and E. Favre, “Optimization of multistage membrane gas separation processes. Example of application to CO2 capture from blast furnace gas,” Journal of Membrane Science, vol. 566, pp. 346 – 366, 2018.
  * [29] A. M. Arias, M. C. Mussati, P. L. Mores, N. J. Scenna, J. A. Caballero, and S. F. Mussati, “Optimization of multi-stage membrane systems for CO2 capture from flue gas,” International Journal of Greenhouse Gas Control, vol. 53, pp. 371 – 390, 2016.
  * [30] V. S. Adi, M. Cook, L. G. Peeva, A. G. Livingston, and B. Chachuat, “Optimization of OSN Membrane Cascades for Separating Organic Mixtures,” in 26th European Symposium on Computer Aided Process Engineering (Z. Kravanja and M. Bogataj, eds.), vol. 38 of Computer Aided Chemical Engineering, pp. 379 – 384, Elsevier, 2016.
  * [31] C. Kunde and A. Kienle, “Global optimization of multistage binary separation networks,” Chemical Engineering and Processing - Process Intensification, vol. 131, pp. 164 – 177, 2018.
  * [32] J. Marriott and E. Sørensen, “The optimal design of membrane systems,” Chemical Engineering Science, vol. 58, no. 22, pp. 4991 – 5004, 2003.
  * [33] J. Marriott and E. Sørensen, “A general approach to modelling membrane modules,” Chemical Engineering Science, vol. 58, no. 22, pp. 4975 – 4990, 2003.
  * [34] S. Weller and W. A. Steiner, “Separation of gases by fractional permeation through membranes,” Journal of Applied Physics, vol. 21, no. 4, pp. 279–283, 1950.
  * [35] S. Weller and W. A. Steiner, “Erratum: Separation of gases by fractional permeation through membranes,” Journal of Applied Physics, vol. 21, no. 12, pp. 1340–1340, 1950.
  * [36] C. Y. Pan, “Gas separation by permeators with high-flux asymmetric membranes,” AIChE Journal, vol. 29, no. 4, pp. 545–552, 1983.
  * [37] M. Ahsan and A. Hussain, “Mathematical modelling of membrane gas separation using the finite difference method,” Pacific Science Review A: Natural Science and Engineering, vol. 18, no. 1, pp. 47 – 52, 2016.
  * [38] X. Feng, J. Ivory, and V. S. V. Rajan, “Air separation by integrally asymmetric hollow-fiber membranes,” AIChE Journal, vol. 45, no. 10, pp. 2142–2152, 1999.
  * [39] R. Qi and M. A. Henson, “Approximate modeling of spiral-wound gas permeators,” Journal of Membrane Science, vol. 121, no. 1, pp. 11 – 24, 1996.
  * [40] R. Qi and M. A. Henson, “Modeling of spiral-wound permeators for multicomponent gas separations,” Industrial & Engineering Chemistry Research, vol. 36, no. 6, pp. 2320–2331, 1997.
  * [41] Y. Kawajir, C. Laird, and A. Wachter, Introduction to Ipopt: A tutorial for downloading, installing, and using Ipopt.  July 2016.
  * [42] J. Viswanathan and I. Grossmann, “A combined penalty function and outer-approximation method for MINLP optimization,” Computers & Chemical Engineering, vol. 14, no. 7, pp. 769 – 782, 1990.
  * [43] R. E. Rosenthal, GAMS – A User’s Guide.  GAMS Development Corporation, Washington, DC, USA, 2014.
  * [44] The Mathworks, Inc., Natick, Massachusetts, MATLAB.
  * [45] J. Xu and R. Agrawal, “Membrane separation process analysis and design strategies based on thermodynamic efficiency of permeation,” Chemical Engineering Science, vol. 51, no. 3, pp. 365 – 385, 1996.
  * [46] J. Wijmans and R. Baker, “The solution-diffusion model: a review,” Journal of Membrane Science, vol. 107, no. 1, pp. 1 – 21, 1995.
  * [47] L. S. White, “Transport properties of a polyimide solvent resistant nanofiltration membrane,” Journal of Membrane Science, vol. 205, no. 1, pp. 191 – 202, 2002.
  * [48] P. Silva, L. G. Peeva, and A. G. Livingston, “Organic solvent nanofiltration (OSN) with spiral-wound membrane elements—Highly rejected solute system,” Journal of Membrane Science, vol. 349, no. 1, pp. 167 – 174, 2010.
  * [49] L. Tranchino, R. Santarossa, F. Carta, C. Fabiani, and L. Bimbi, “Gas separation in a membrane unit: Experimental results and theoretical predictions,” Separation Science and Technology, vol. 24, no. 14, pp. 1207–1226, 1989.
  * [50] H. K. Lonsdale, U. Merten, and M. Tagami, “Phenol transport in cellulose acetate membranes,” Journal of Applied Polymer Science, vol. 11, no. 9, pp. 1807–1820, 1967.
  * [51] D. Paul, “Reformulation of the solution-diffusion theory of reverse osmosis,” Journal of Membrane Science, vol. 241, no. 2, pp. 371 – 386, 2004.
  * [52] H. Zimmermann, Propene.  Wiley-VCH Verlag GmbH & Co. KGaA, 2013.
  * [53] J. R. Alcántara-Avila, F. I. Gómez-Castro, J. G. Segovia-Hernández, K.-I. Sotowa, and T. Horikawa, “Optimal design of cryogenic distillation columns with side heat pumps for the propylene/propane separation,” Chemical Engineering and Processing: Process Intensification, vol. 82, pp. 112 – 122, 2014.
  * [54] F. Fábrega, J. Rossi, and J. d’Angelo, “Exergetic analysis of the refrigeration system in ethylene and propylene production process,” Energy, vol. 35, no. 3, pp. 1224 – 1231, 2010.
  * [55] C. W. Colling, G. A. Huff Jr, and J. V. Bartels, “Processes using solid perm-selective membranes in multiple groups for simultaneous recovery of specified products from a fluid mixture,” Jan. 2004.  US Patent 20040004040 A1.
  * [56] A. Górak and H. Schoenmakers, Distillation Operation and Applications.  Boston: Academic Press, 2014.
  * [57] J. Klabunde, C. Bischoff, and A. J. Papa, Propanols.  Wiley-VCH Verlag GmbH & Co. KGaA, 2018.
  * [58] D. R. Summers, P. J. McGuire, M. R. Resetarits, C. E. Graves, S. E. Harper, and S. J. Angelino, “High-capacity trays debottleneck Texas C3 splitter,” Oil and Gas Journal, vol. 93, no. 45, 1995.
  * [59] M. O. Daramola, Xylenes: Synthesis, Characterization and Physicochemical Properties.  Nova Science Publishers, Incorporated, 2013.
  * [60] C. Perego and P. Pollesel, “Chapter 2 - advances in aromatics processing using zeolite catalysts,” in Advances in Nanoporous Materials (S. Ernst, ed.), vol. 1 of Advances in Nanoporous Materials, pp. 97 – 149, Elsevier, 2010.
  * [61] M. T. Ashraf, Analysis and Optimization of p-Xylene Production Process.  PhD thesis, 2013.
  * [62] T.-C. Tsai, S.-B. Liu, and I. Wang, “Disproportionation and transalkylation of alkylbenzenes over zeolite catalysts,” Applied Catalysis A: General, vol. 181, no. 2, pp. 355 – 398, 1999.
  * [63] D.-Y. Koh, B. A. McCool, H. W. Deckman, and R. P. Lively, “Reverse osmosis molecular differentiation of organic liquids using carbon molecular sieve membranes,” Science, vol. 353, no. 6301, pp. 804–807, 2016.
  * [64] R. Chafin, J. S. Lee, and W. J. Koros, “Effects of casting and post casting annealing on xylene isomer transport properties of Torlon® 4000T films,” Polymer, vol. 51, no. 15, pp. 3462 – 3471, 2010.
  * [65] L. M. Robeson, “The upper bound revisited,” Journal of Membrane Science, vol. 320, no. 1, pp. 390 – 400, 2008.