Robust Multi-Objective Optimization of Chromatographic Rare Earth Element Separation

Show more

Robust Optimization

1. Introduction

Rare earth elements (REE) are extensively used in modern technological industries [1] - [7] , and are considered as strategic commodities in many countries [4] [6] [7] [8] . REE minerals with varying compositions are found at deposits throughout the world [1] [2] [3] [4] [5] [7] [9] , though most of the global REE supply comes from only a few sources [1] . By using liquid-liquid extraction methods, the elements are separated from the minerals and upgraded to suitable purity levels for commercial applications [1] [2] [3] [4] [7] . However, several experimental and model-based studies have accentuated chromatography as an alternative method with considerable benefits [1] [5] [10] - [17] . This study is intended as a contribution to the work of developing chromatography as an REE separation method, and focuses on preparative chromatographic batch separation of the middle REE group, samarium (Sm), europium (Eu) and gadolinium (Gd), where Eu is the target product due to its’ higher commercial value. It is a continuation of an experimental optimization study [11] , which was followed by a model based multi-objective optimization (MOO) study [18] where a process optimization strategy was presented. The current work complements the previous studies by introducing robust multi-objective optimization. This is utilized for mapping the required operation point changes, needed to keep the number of failed batches at an acceptable level when a certain level of process disturbance is introduced, as well as evaluating the performance change that is accounted for when formulating a robust counterpart problem.

Since mathematical modelling offers a cost efficient and powerful approach for assessing preparative chromatography [10] [13] [14] [19] - [25] , it has been employed as the preferred tool for evaluating and optimizing the chromatographic system in this study. The optimization of a chromatography system is ordinarily cast in a bi-level framework [25] : 1) the upper level that administer the effects of the decision variables, such as load and elution gradient, that governs the chromatogram, and 2) the lower level that establishes the pooling strategy for deciding the product pooling cut-times. A MOO method is needed when competing optimization objectives, such as productivity and yield, are considered. The MOO strategy in this work, as presented in [18] , provides with an intact optimization objective for all levels of the bi-level optimization and firm objective values when evaluating the multi-objective optimization problem (MOP). However, the nominal solution for a MOP is not necessarily robust, and even small process disturbances may cause process failure. This calls for a transformation of the MOP into its’ robust counterpart problem [26] [27] .

Approaches for introducing robustness to process optimizations are readily available in literature, and robust optimization for chromatography in particular are, amongst others, presented in [20] [24] [27] - [32] . The preferred robust optimization method used in this work is similar to the methods presented in [29] [33] , with the difference of that this study exclusively utilizes a stochastic method to obtain the model responses of the introduced process disturbances.

The main focus of this study was to perform a robust multi-objective optimization of chromatographic REE separation. In this context, an experimentally validated process model from a previous study [18] was used to generate the process system response. The results from the robust multi-objective optimization were used to assess the robustness of the system, chart the process parameter changes when robustness was introduced, and evaluate the expected loss of performance when robustness was applied.

2. Chromatography Model

2.1. Process Description

The chromatography model in this work is based on an experimental study of batch chromatography separation of Sm, Eu and Gd [11] . The experimental study utilized an Agilent 1200 series HPLC system (Agilent Technologies, Waldbronn, Germany) and a Kromasil M3 (Eka, Bohus, Sweden) column with the dimensions 150 × 4.6 mm. An inductively coupled plasma mass spectrometry (ICP-MS) system (Agilent Technologies, Tokyo, Japan) was used for in-line post column REE detection. The column stationary phase was made of spherical, C18 coated, silica particles with a diameter of 16 μm and a pore size of 100 Å. Each column had a ligand concentration of 342 mM Bis (2-ethylhexyl) phosphoric acid (Sigma-Aldrich, St. Louis, USA), and the elution gradient concentration was set to vary between 6% - 13% (vol) of 7 M nitric acid over a gradient length of 5 column volumes (CV).

2.2. Chromatography Model

The chromatography model in this study has been presented in a previous publication [18] , and is reproduced here for the purpose of clarity. The column separation was described through a kinetic dispersive model [34] with a Langmuir mobile phase modulator isotherm [13] [14] [20] [21] . The model equations, defined in the spatial, $z\in \left[{z}_{0}\mathrm{,}{z}_{f}\right]$ , and temporal, $t\in \left[{t}_{0}\mathrm{,}{t}_{f}\right]$ , domains are given by:

$\frac{\partial {c}_{\alpha}}{\partial t}=-\frac{\partial}{\partial z}\left({c}_{\alpha}{v}_{\text{int}}-{D}_{\text{app,}\alpha}\frac{\partial {c}_{\alpha}}{\partial z}\right)-\frac{\left(1-{\epsilon}_{c}\right)}{{\epsilon}_{c}+\left(1-{\epsilon}_{c}\right){\epsilon}_{p}}\frac{\partial {q}_{\alpha}}{\partial t}\mathrm{,}$ (1)

$\frac{\partial {q}_{\alpha}}{\partial t}={k}_{\text{kin}\mathrm{,}\alpha}\left({c}_{\alpha}{K}_{eq\mathrm{,}\alpha}{q}_{\mathrm{max},\alpha}\left[1-{\displaystyle \underset{\gamma \in \left\{\text{Sm,Eu,Gd}\right\}}{\sum}}\frac{{q}_{\gamma}}{{q}_{\text{max}\mathrm{,}\gamma}}\right]-{q}_{\alpha}{c}_{S}^{{\nu}_{\alpha}}\right)\mathrm{,}$ (2)

where ${c}_{\alpha}$ and ${q}_{\alpha}$ are the mobile and solid phase concentration of component $\alpha \in \left\{\text{Sm},\text{Eu},\text{Gd},\text{S}\right\}$ , ${v}_{\text{int}}$ is the quotient of superficial velocity over total porosity, ${D}_{\text{app,\alpha}}$ the apparent dispersion coefficient, and ${\epsilon}_{c}$ and ${\epsilon}_{p}$ the column and particle void fractions. Here, ${c}_{S}$ denotes the concentration of the modifier (i.e. nitric acid), ${k}_{\text{kin}\mathrm{,}\alpha}$ a parameter describing the kinetics, ${K}_{eq\mathrm{,}\alpha}$ the equilibirum constant regarding adsorption and desorption, ${\nu}_{\alpha}$ a parameter describing the ion-exchange characteristics, and ${q}_{\text{max}\mathrm{,}\alpha}$ the maximum concentration of adsorbed components. The model does not consider modifier ions on the solid phase, therefore Equation (2) and its associated part in Equation (1) are omitted (i.e. $\partial {q}_{\alpha}/\partial t\equiv 0$ ) when $\alpha =S$ . Equation (1) is complemented with Danckwert boundary conditions [35] :

$\begin{array}{l}{c}_{\alpha}\left(t\mathrm{,}{z}_{0}\right){v}_{\text{int}}-{D}_{\text{app,\alpha}}\frac{\partial {c}_{\alpha}}{\partial z}\left(t\mathrm{,}{z}_{0}\right)\\ =(\begin{array}{ll}{c}_{\text{load}\mathrm{,}\alpha}{v}_{\text{int}}\Pi \left(t\mathrm{,}{t}_{0}\mathrm{,}\Delta {t}_{\text{load}}\right)\hfill & \text{\hspace{0.05em}}\text{if}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\alpha \in \left\{\text{Sm,Eu,Gd}\right\}\mathrm{,}\hfill \\ {c}_{\text{mix}\mathrm{,}S}{v}_{\text{int}}\hfill & \text{\hspace{0.05em}}\text{if}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\alpha =S\mathrm{,}\hfill \end{array}\end{array}$ (3)

$\frac{\partial {c}_{\alpha}}{\partial z}\left(t\mathrm{,}{z}_{f}\right)=\mathrm{0,}\text{\hspace{1em}}\forall \alpha \in \left\{\text{Sm},\text{Eu},\text{Gd},\text{S}\right\}$ (4)

where ${c}_{\text{load}\mathrm{,}\alpha}$ is the injected load concentration, and $\Pi \left(t\mathrm{,}{t}_{0}\mathrm{,}\Delta {t}_{\text{load}}\right)\in \left\{\mathrm{0,1}\right\}$ a rectangular function in the temporal horizon $\left[{t}_{0}\mathrm{,}\Delta {t}_{\text{load}}\right]$ . The dynamics of the modifier concentration in the upstream mixing tank, ${c}_{\text{mix}\mathrm{,}S}$ , are given by:

$\frac{\text{d}{c}_{\text{mix}\mathrm{,}S}}{\text{d}t}=\frac{1}{{\tau}_{\text{mix}}}\left(u\left(t\right)-{c}_{\text{mix}\mathrm{,}S}\right)\mathrm{,}$ (5)

$u\left(t\right)=(\begin{array}{ll}{u}_{0}\mathrm{,}\hfill & \text{\hspace{0.05em}}\text{if}\text{\hspace{0.17em}}\text{\hspace{0.05em}}t\le \Delta {t}_{\text{load}}+\Delta {t}_{\text{wash}}\mathrm{,}\hfill \\ {u}_{0}+t-\left(\Delta {t}_{\text{load}}+\Delta {t}_{\text{wash}}\right)\mathrm{,}\hfill & \text{\hspace{0.05em}}\text{if}\text{\hspace{0.17em}}\text{\hspace{0.05em}}t>\Delta {t}_{\text{load}}+\Delta {t}_{\text{wash}}\mathrm{,}\hfill \end{array}$ (6)

where ${\tau}_{\text{mix}}$ is the residence time, $u$ is the elution gradient described by the initial value, ${u}_{0}$ , and the slope of the linear elution gradient, $\Delta u$ , expressed as

$\Delta u=\frac{{u}_{f}-{u}_{0}}{{t}_{f}-\left(\Delta {t}_{\text{load}}+\Delta {t}_{\text{wash}}\right)}$ .

The first-order spatial derivative in Equation (1) was approximated using a method-of-lines and finite volume method with 100 grid points where ${z}_{k}=k\Delta z$ is the discretized spatial coordinate and $k\in \left[1,100\right]$ . The first order derivative was approximated as a two-point backward difference, and the second-order derivative was approximated as a three-point central difference. The model parameter values from [18] were used in this study, and are given in Table 1.

3. Robust Multi-Objective Optimization

3.1. Multi-Objective Optimization

The multi-objecitve optimization problem formulation in this work resembles that of a previous study [18] , with the difference of that this work only considers the two competing objectives yield and productivity. This is in order to benchmark the proposed robust optimization method with a previous bi-objective robust optimization method [30] . In this study, the column outlet concentration profile, ${c}_{\alpha}\left(t\mathrm{,}{z}_{f}\right)$ where $\alpha \in \left\{\text{Sm},\text{Eu},\text{Gd}\right\}$ , was used for evaluation of the competing objective functions yield, ${Y}_{\alpha}$ and productivity, ${P}_{\alpha}$ for the target component Eu. The objective functions for the collected component, $\alpha $ , between the

Table 1. Model parameter values.

cut-times $\left[{t}_{c}\mathrm{,}{t}_{f}\right]$ are defined as:

${\delta}_{\text{load}\mathrm{,}\alpha}\frac{\text{d}{Y}_{\alpha}}{\text{d}t}={c}_{\alpha}\left(t\mathrm{,}{z}_{f}\right){v}_{\text{int}}{A}_{c}\Pi \left(t\mathrm{,}{t}_{c}\mathrm{,}{t}_{f}\right)\mathrm{,}$ (7)

$\frac{\text{d}{P}_{\alpha}}{\text{d}t}=\frac{1}{{V}_{c}}\frac{1}{{t}_{f}+{t}_{r}}{\delta}_{\text{load}\mathrm{,}\alpha}\frac{\text{d}{Y}_{\alpha}}{\text{d}t}\mathrm{,}$ (8)

where
${\delta}_{\text{load}\mathrm{,}\alpha}={c}_{\text{load}\mathrm{,}\alpha}{v}_{\text{int}}{A}_{c}\Delta {t}_{\text{load}}$ is the total amount of injected sample,
${A}_{c}$ and V_{c} the column cross-sectional area and volume, and
${t}_{r}=2{V}_{c}{\stackrel{\dot{}}{Q}}^{-1}$ the regeneration and re-equilibration time following the final cut-time. Thus, the objective becomes to determine an optimal elution gradient,
$u$ , batch load,
${\delta}_{\text{load}\mathrm{,}\alpha}$ , and pooling cut-times,
$\left[{t}_{c}\mathrm{,}{t}_{f}\right]$ , that maximizes
${Y}_{\alpha}\left({t}_{f}\right)$ and
${P}_{\alpha}\left({t}_{f}\right)$ , while fulfilling the target component purity constraint given by:

${X}_{\alpha}\left({t}_{f}\right)=\frac{{\delta}_{\text{load}\mathrm{,}\alpha}{Y}_{\alpha}\left({t}_{f}\right)}{{\displaystyle \underset{\beta \in \left\{\text{Sm},\text{Eu},\text{Gd}\right\}}{\sum}}{\delta}_{\text{load}\mathrm{,}\beta}{Y}_{\beta}\left({t}_{f}\right)}\mathrm{,}$ (9)

where the numerator is the captured amount of the target component in $\left[{t}_{c}\mathrm{,}{t}_{f}\right]$ and the denominator represents the total amount of captured components.

The weighted sum scalarization method was used to combine the objectives in Equations (7)-(8) to a single objective with the weight for productivity defined as $\omega \in \left[\mathrm{0,1}\right]$ , and the weight for yield defined as $1-\omega $ . The decision variables are the free operating parameters, i.e. $\Delta {t}_{\text{load}}\mathrm{,}{t}_{c}\mathrm{,}{t}_{f}\mathrm{,}{u}_{0}$ and ${u}_{f}$ , which in turn determine the trajectories $x=\left({c}_{\alpha}\mathrm{(}t\mathrm{,}{z}_{k}\mathrm{),}{c}_{S}\left(t\mathrm{,}{z}_{k}\right)\mathrm{,}{c}_{\text{mix}\mathrm{,}S}\left(t\right)\mathrm{,}{q}_{\alpha}\left(t\mathrm{,}{z}_{k}\right)\mathrm{,}{P}_{\text{Eu}}\left(t\right)\mathrm{,}{Y}_{\text{Eu}}\left(t\right)\mathrm{,}{X}_{\text{Eu}}\left(t\right)\right)$ . The resulting optimization problem can then be set in the framework for min-min optimal control:

$\text{min}\text{.}\text{\hspace{1em}}\text{\hspace{0.05em}}-\left(\omega {\displaystyle \underset{{t}_{0}}{\overset{{t}_{f}}{\int}}}\frac{\text{d}{P}_{\text{Eu}}}{\text{d}t}\text{d}t+\left(1-\omega \right){\displaystyle \underset{{t}_{0}}{\overset{{t}_{f}}{\int}}}\frac{\text{d}{Y}_{\text{Eu}}}{\text{d}t}\text{d}t\right)\mathrm{,}$ (10a)

$\text{w}\text{.r}\text{.t}\text{.}\text{\hspace{1em}}p=\left(\Delta {t}_{\text{load}}\mathrm{,}{u}_{0}\mathrm{,}{u}_{f}\right)\in {\mathbb{R}}^{3}\mathrm{,}$

$\text{s}\text{.t}\text{.}\text{\hspace{1em}}\text{\hspace{1em}}{p}_{L}\le p\le {p}_{U}\mathrm{,}$ (10b)

$\left(x\mathrm{,}{t}_{c}\mathrm{,}{t}_{f}\right)=\text{argmin}\text{.}-\left(\omega {\displaystyle \underset{{t}_{0}}{\overset{{t}_{f}}{\int}}}\frac{\text{d}{P}_{\text{Eu}}}{\text{d}t}\text{d}t+\left(1-\omega \right){\displaystyle \underset{{t}_{0}}{\overset{{t}_{f}}{\int}}}\frac{\text{d}{Y}_{\text{Eu}}}{\text{d}t}\text{d}t\right)\mathrm{,}\text{\hspace{0.05em}}$ (10c)

$\text{w}\text{.r}\text{.t}\text{.}\text{\hspace{1em}}\left({t}_{c}\mathrm{,}{t}_{f}\right)\in {\mathbb{R}}^{2}\mathrm{,}$

$\text{s}\text{.t}\text{.}\text{\hspace{1em}}\stackrel{\dot{}}{x}=F\left(t\mathrm{,}x\left(t\right)\mathrm{,}{t}_{c}\mathrm{,}{t}_{f}\mathrm{,}p\right)\mathrm{,}\text{\hspace{1em}}x\left({t}_{0}\right)={x}_{0}\mathrm{,}$ (10d)

${X}_{\text{Eu}\mathrm{,}L}-{X}_{\text{Eu}}\left({t}_{f}\right)\le \mathrm{0,}$ (10e)

${t}_{c\mathrm{,}L}\le {t}_{c}\le {t}_{c\mathrm{,}U}\mathrm{,}\text{\hspace{1em}}{t}_{f\mathrm{,}L}\le {t}_{f}\le {t}_{f\mathrm{,}U}\mathrm{,}$ (10f)

$\forall t\in \left[{t}_{0}\mathrm{,}{t}_{f}\right]\mathrm{,}\text{\hspace{0.05em}}\text{\hspace{1em}}\forall z\in \left[{z}_{0}\mathrm{,}{z}_{f}\right]\mathrm{.}$

The solution of Equation (10) will result in a Pareto front situated on the boundary of the feasible region. This implies that a process disturbance, however slight, can cause a batch failure in terms of not meeting the purity constraint [26] [29] . In order to account for such disturbances it is necessary to formulate a robust counterpart problem, which in this study will be accomplished by introducing a back-off term to the purity inequality constraint in Equation (10e).

3.2. Robust Multi-Objective Optimization Problem Formulation

In order to formulate a robust counterpart of Equation (10), we consider a set of bounded distributed disturbances, $\stackrel{\u02dc}{p}$ , on the free operating parameters, $p$ (i.e $\Delta {t}_{\text{load}}\mathrm{,}{u}_{0}$ and ${u}_{f}$ ), and define ${\stackrel{\u02dc}{X}}_{\text{Eu}}$ as the cumulative purity distribution of the model responses that are produced from $\stackrel{\u02dc}{p}$ . A purity constraint back-off term, ${X}_{\text{BF}}$ , is introduced in order to make the purity constraint robust with respect to the disturbances. The back-off term can essentially be seen as a safety margin that amplifies the purity inequality constraint in Equation (10e) so that the purity requirement, ${X}_{\text{Eu,}L}$ , still can be met for the considered set of bounded disturbances. The success rate is defined as the fraction of batches in the disturbance set that fulfil the purity requirement, ${X}_{\text{Eu,}L}$ , and ${\Phi}_{{X}_{\text{Eu}}}$ signifies the desired success rate. The following robust counterpart of Equation (10) is then given by:

$\text{min}\text{.}\text{\hspace{1em}}\text{\hspace{0.05em}}{\displaystyle \underset{-\infty}{\overset{{X}_{\text{Eu}\mathrm{,}L}+{X}_{\text{BF}}}{\int}}}{\stackrel{\u02dc}{X}}_{\text{Eu}}\text{d}{X}_{\text{Eu}}-{\Phi}_{{X}_{\text{Eu}}}\mathrm{,}$ (11a)

$\text{w}\text{.r}\text{.t}\text{.}\text{\hspace{1em}}{X}_{\text{BF}}\mathrm{,}$

$\text{s}\text{.t}\text{.}\text{\hspace{1em}}\stackrel{\dot{}}{\stackrel{\u02dc}{x}}=\stackrel{\u02dc}{F}\left(t\mathrm{,}x\left(t\right)\mathrm{,}{t}_{c}\mathrm{,}{t}_{f}\mathrm{,}\stackrel{\u02dc}{p}\right)\mathrm{,}$ (11b)

$\stackrel{\u02dc}{p}~N\left(p\mathrm{,}{\sigma}_{p}^{2}\right)\mathrm{,}$ (11c)

$p=\text{argmin}\text{.}-\left(\omega {\displaystyle \underset{{t}_{0}}{\overset{{t}_{f}}{\int}}}\frac{\text{d}{P}_{\text{Eu}}}{\text{d}t}\text{d}t+\left(1-\omega \right){\displaystyle \underset{{t}_{0}}{\overset{{t}_{f}}{\int}}}\frac{\text{d}{Y}_{\text{Eu}}}{\text{d}t}\text{d}t\right)\mathrm{,}$ (11d)

$\text{w}\text{.r}\text{.t}\text{.}\text{\hspace{1em}}p=\left(\Delta {t}_{\text{load}}\mathrm{,}{u}_{0}\mathrm{,}{u}_{f}\right)\in {\mathbb{R}}^{3}\mathrm{,}$

$\text{s}\text{.t}\text{.}\text{\hspace{1em}}{p}_{L}\le p\le {p}_{U}\mathrm{,}$ (11e)

$\left(x\mathrm{,}{t}_{c}\mathrm{,}{t}_{f}\right)=\text{argmin}\text{.}-\left(\omega {\displaystyle \underset{{t}_{0}}{\overset{{t}_{f}}{\int}}}\frac{\text{d}{P}_{\text{Eu}}}{\text{d}t}\text{d}t+\left(1-\omega \right){\displaystyle \underset{{t}_{0}}{\overset{{t}_{f}}{\int}}}\frac{\text{d}{Y}_{\text{Eu}}}{\text{d}t}\text{d}t\right)\mathrm{,}$ (11f)

$\text{w}\text{.r}\text{.t}\text{.}\text{\hspace{1em}}\left({t}_{c}\mathrm{,}{t}_{f}\right)\in {\mathbb{R}}^{2}\mathrm{,}$

$\text{s}\text{.t}\text{.}\text{\hspace{1em}}\stackrel{\dot{}}{x}=F\left(t\mathrm{,}x\left(t\right)\mathrm{,}{t}_{c}\mathrm{,}{t}_{f}\mathrm{,}p\right)\mathrm{,}\text{\hspace{1em}}x\left({t}_{0}\right)={x}_{0}\mathrm{,}$ (11g)

$\left({X}_{\text{Eu}\mathrm{,}L}+{X}_{\text{BF}}\right)-{X}_{\text{Eu}}\left({t}_{f}\right)\le \mathrm{0,}$ (11h)

${t}_{c\mathrm{,}L}\le {t}_{c}\le {t}_{c\mathrm{,}U}\mathrm{,}\text{\hspace{1em}}{t}_{f\mathrm{,}L}\le {t}_{f}\le {t}_{f\mathrm{,}U}\mathrm{,}$ (11i)

$\forall t\in \left[{t}_{0}\mathrm{,}{t}_{f}\right]\mathrm{,}\text{\hspace{1em}}\forall z\in \left[{z}_{0}\mathrm{,}{z}_{f}\right]\mathrm{.}$

A decomposition strategy was adopted to transform the robust MOP into three levels: 1) the upper-level optimization problem given by Equations (11a)-(11c) with respect to ${X}_{\text{BF}}$ , 2) the mid-level optimization problem given by Equations (11d)-(11e) with respect to p, and 3) the lower-level optimization problem given by Equations (11f)-(11i) and constrained by the ODE system, $F$ , governed by Equations (1), (2), (5), (7), (8), (9). Essentially, Equation (11) was solved by using the simulated system response, $\stackrel{\u02dc}{x}$ , for an uncertainty set of the free operating parameters, $\stackrel{\u02dc}{p}$ , to evaluate the cumulative distribution function of ${\stackrel{\u02dc}{X}}_{\text{Eu}}$ . The back-off term, ${X}_{\text{BF}}$ , in the purity inequality constraint, Equation (11h), was then incrementally increased to gain more successful batches in ${\stackrel{\u02dc}{X}}_{\text{Eu}}$ , and thereby achieving a more robust process. This procedure was repeated iteratively until Equation (11a) was fulfilled, at which point Equation (11) was considered to be solved.

In this study, the desired success rate, ${\Phi}_{{X}_{\text{Eu}}}$ , was set to 0.95, and the decision variable boundaries are presented in Table 2.

3.3. Optimization Method

The robust optimization method in this study is similar to the methods presented in [29] [33] . The main difference is that the model responses of the introduced disturbances in this study are obtained stochastically, instead of alternatively utilizing a deterministic approach through linearization of the uncertainty set. The stochastic approach has the benefit of being more straightforward, at the expense of an increased demand of computation power. This increased demand was accommodated for by using a parallel computing

Table 2. Decision variables.

methodology as described in [19] .

As a first step, the nominal and non-robust Pareto front was obtained by solving the MOP as defined in Equation (10). This was carried out through MATLAB’s fmincon function with a sequential quadratic programming algorithm, the BFGS formula for updating the approximation of the Hessian matrix, and central differences to estimate the gradient of the objective function and constraint functions. Then an uncertainty set, $\stackrel{\u02dc}{p}$ , with a normal distribution, assuming no covariance between the free operating parameters $p$ , a standard deviation $\sigma $ , and sampling size of 10.000 was obtained via MATLAB’s lhsnorm function. The uncertainty set was applied to the investigated operating points on the nominal Pareto front, and the model responses were used to evaluate the cumulative purity distribution, ${\stackrel{\u02dc}{X}}_{\text{Eu}}$ , of the uncertainty set.

Then, an initial investigation of the back-off terms impact on ${\stackrel{\u02dc}{X}}_{\text{Eu}}$ was conducted by creating new Pareto fronts with an incrementally increased back-off and observing how ${\stackrel{\u02dc}{X}}_{\text{Eu}}$ changes when $\stackrel{\u02dc}{p}$ is applied to the investigated points on the new Pareto fronts. At this stage, it is of particular interest to investigate how the fraction of batches that fulfil the purity requirement in the perturbed set, changes with an increased back-off. This provides an estimate of the required back-off to meet a certain success rate for a given purity constraint.

The required back-off for a given point on the nominal Pareto front was obtained by applying MATLAB’s fminbnd function on the upper level of the robust counterpart problem in Equation (11), with suitable boundaries obtained from the previous back-off investigation. The mid- and lower-level optimization problems in Equation (11) were solved by MATLAB’s fmincon function with a sequential quadratic programming algorithm, the BFGS formula for updating the approximation of the Hessian matrix, and central differences to estimate the gradient of the objective function and constraints. The procedure comprises an evaluation of the cumulative distribution function of ${\stackrel{\u02dc}{X}}_{\text{Eu}}$ based on $\stackrel{\u02dc}{x}$ and $\stackrel{\u02dc}{p}$ , as obtained from the mid- and lower-level optimization problem for a given initial ${X}_{\text{BF}}$ . ${X}_{\text{BF}}$ is then varied for the upper level optimization problem through MATLAB's fminbnd function, resulting in new $\stackrel{\u02dc}{x}$ , $\stackrel{\u02dc}{p}$ and cumulative distribution functions of ${\stackrel{\u02dc}{X}}_{\text{Eu}}$ to be evaluated. This continues until a ${X}_{\text{BF}}$ that produces a cumulative distribution function of ${\stackrel{\u02dc}{X}}_{\text{Eu}}$ corresponding to the desired success rate ${\Phi}_{{X}_{\text{Eu}}}$ is obtained.

Optimization Method Benchmarking

The proposed optimization method was compared with a previous optimization method [30] that focuses on the nominal Pareto front and optimizes the pooling time horizon for each investigated point on the front so that the purity requirement is met for a given uncertainty set. The main difference between the two methods is that the proposed method will find new optimal operation points by changing the free operating parameters, and achieve robustness by increasing the purity requirement back-off for each point on the Pareto front, whereas the previous method keeps the decision variables from the nominal Pareto front intact, with the exception of the cut-times that are optimized to find a fixed pooling time horizon that will fulfil the purity requirement for the entire uncertainty set.

3.4. Robust Multiobjective Optimization Results

The robust multiobjective optimizations of the studied system were carried out with a product purity requirement, ${X}_{\text{Eu},L}$ , of 0.95 and 0.99 respectively, and the target success rate, ${\Phi}_{{X}_{\text{Eu}}}$ , was set to 0.95. The perturbed process parameters were the injected load concentration, ${c}_{\text{load}\mathrm{,}\alpha}$ , and the modifier concentration in the upstream mixing tank, ${c}_{\text{mix}\mathrm{,}S}$ . Several values for the uncertainty set standard deviation, $\sigma $ , were investigated. However, a standard deviation exceeding 0.01 did not result in achieving the target success rate even when the back-off was set to the maximum limit, i.e. ${X}_{\text{Eu}\mathrm{,}L}+{X}_{\text{BF}}=1$ . Therefore, only results for $\sigma =0.01$ are presented in this study. This should be interpreted as that the system is very un-robust, and sensitive to even the slightest process disturbances. However, this was expected since the studied elements are extremely similar in both chemical and physical properties, resulting in a minute separation selectivity which in turn makes the separation very difficult and unforgiving towards process perturbations.

The nominal un-robust Pareto fronts are presented by the outermost fronts in Figure 1, and it is shown how the Pareto front decreases with an increased back-off on the purity requirement. It should be noted that we only present the impact of an increased back-off on four Pareto points to avoid overtly crowded figures. The objective weight, $\omega $ , for these points correspond to 1 (i.e. productivity as single objective), 0.5, 0.3 and 0 (i.e. yield as single objective) respectively.

Figure 1. The nominal Pareto fronts are presented by the solid lines for a purity requirement of 0.95 in (a) and 0.99 in (b). The cross marks indicate how a Pareto point, with the weight $\omega $ , changes with an increased back-off. The dashed lines indicate the Pareto front outlines for an increasing back-off, and it can be seen that the Pareto front decreases as the back-off is increased.

The results from the initial investigation on how the success rate, ${\Phi}_{{X}_{\text{Eu}}}$ , for the investigated points on the nominal Pareto front increases with an increased back-off are shown in Figure 2. The figure provides with an estimation of the required back-off to achieve the desired success rate for a given disturbance set, and it can be seen that an objective leaning more towards yield (i.e. $\omega $ decreases) results in a lower success rate for a given back-off.

It should be noted that the Pareto point corresponding to yield as a single objective, i.e. $\omega =0$ , has been omitted since the point for maximum yield is considered as an undesirable operating point due to the drastic decrease in productivity.

It is somewhat counter intuitive that the success rate should decrease with an increased objective weight for yield, since a higher yield typically is associated with an increased peak separation which in turn should result in an increased robustness. The decrease of robustness can be explained by observing how the decision variables change with an increased back-off for the 0.95 purity requirement case in Figure 3, where the pooling cut-time trends become very interesting. The decision variable trends show that the initial elution concentration and elution gradient slope are quite similar as long as productivity is a part of the weighted objective. However, a higher productivity is favoured by a larger batch load, a pooling horizon occurring earlier in the chromatogram (i.e. first and last pooling cuts occur earlier) and a smaller pooling volume. The increased batch load is reasonable since it will allow for a higher productivity due to an increased throughput. The early first cut comes from that a higher batch load will capacitate the elements to start eluting earlier. The earlier final cut makes the cycle time shorter, which is favourable for productivity, but it is also a trade of in terms of decreased yield.

Figure 2. Results from the investigation of how the success rate increases with an increased back-off for a purity requirement of 0.95 in (a) and 0.99 in (b). The dashed line indicates the target success rate, ${\Phi}_{{X}_{\text{Eu}}}$ , and helps to provide an initial estimation of the required back-off, ${X}_{\text{BF}}$ , for a Pareto point with the objective weight $\omega $ .

Figure 3. Plots of decision variable changes due to an increasing back-off for Pareto points with different objective weights, $\omega $ , and a purity requirement of 0.95. (a) Batch load, (b) Initial elution concentration, (c) Elution gradient slope, (d) First cut time, (e) Final cut time, (f) Pooling volume.

This has the important implication of that a high objective weight on productivity will result in pooling cut times occurring closer to the Eu elution peak centre and farther away from the neighbouring peaks. When a higher yield is desired, the pooling horizon will increase in order to capture more of the target molecules, and this will move the pooling close to, and even into, the neighbouring elution peaks as long as the purity requirement is met. For this reason, a higher weight on yield will demand a higher back-off on purity in order to meet the desired success rate. This is due to that when a perturbation is introduced, the neighbouring peaks may move closer to, and even intrude, the pooling horizon, and a higher purity requirement will move the pooling cut times farther away from the neighbouring peaks. The farther away the cut times are from the neighbouring peaks in the nominal case, the higher disturbance can be tolerated since there is more room available for the neighbouring peaks to move before they impact the purity of the target peak. This is illustrated in Figure 4 where we can see a case with high productivity (smaller pooling horizon) and a case with high yield (larger pooling horizon) and observe how the introduced process disturbances make the neighbouring peaks creep into the pooling horizon to a larger extent for the high yield case.

The results from the robust optimizations can be seen in Figure 5 where the nominal un-robust Pareto front is plotted along with the robust Pareto front produced by the proposed method, and the front from the benchmark method. It should be noted that both methods provide with robust operating points that

Figure 4. Sections of chromatograms with focus on the collection of the Eu product pool. The purity requirement, ${X}_{\text{Eu}}$ , was set to 0.95 for the productivity objective weights $\omega =1$ (a) and $\omega =0.3$ (b). The pooling cut times are indicated by the dotted lines, and the elution order is Sm, Eu and Gd. The unit for the nitric acid elution gradient (black dashed line) on the vertical axis is mol/l. The shaded areas indicate the span of concentration profile variations due to process disturbances with $\sigma =0.01$ , and the solid black lines indicate the concentration profiles for the nominal case. The chromatograms demonstrate that an operation point with a higher objective weight for productivity (a), is more robust than an operation point with a lower weight (b). This can be seen by observing how the larger pooling horizon in (b) allows for more collection of the neighbouring elements when disturbances are introduced, and thereby causing an increased number of batches with failed purity requirement. This is particularly noticeable for the Gd peak which intrudes the collected pool to a larger extent when process disturbances are introduced in (b).

Figure 5. Pareto fronts resulting from the optimizations with a purity requirement of 0.95 (left) and 0.99 (right). (a) indicates the nominal Pareto front, (b) the robust Pareto front according to the proposed method and (c) the robust front from the benchmark method. The dots indicate the system response of the distributed uncertainty set associated with the respective Pareto points. The dashed lines indicate the loss of productivity for a given yield when robustifying the nominal Pareto front according to the proposed method.

handle the given process disturbances satisfactorily. However, the proposed method should be favoured since it produces a robust Pareto front with higher objective values compared to the benchmark method, which implies that the benchmark method should be considered more restrictive. Further, the benchmark method generates operating points that cannot be considered Pareto optimal, which is the case for the points with $\omega =1$ on front (c) in Figure 5. For the sake of fairness and as mentioned in [36] , these points should be disregarded when generating a robust Pareto front with the benchmark method.

Applying robustness to a point on the nominal Pareto front with a given $\omega $ will result in a change of both productivity and yield, as can be seen in Figure 1, and this makes the evaluation of performance loss when introducing robustness slightly ambiguous. In order to resolve this, the productivity for a given yield on the nominal Pareto front is compared to the productivity on the robust Pareto front given the same yield, as indicated by red dashed lines in Figure 5. The loss of productivity is presented in Table 3, and it shows that a productivity loss in the range of 10% - 20% can be expected. It should be mentioned that the robustness can be increased further during operation by applying a variable pooling cut time control strategy as described in [31] , and thereby decrease the loss of productivity.

The required back-off, that meets the robustness requisite according to the proposed method, is presented in Table 4 for the investigated Pareto points. The results show that a higher purity requirement calls for a larger back-off term to achieve a robust Pareto point, and a lower objective weight on productivity also necessitates an increased back-off.

Table 3. Loss of productivity.

Table 4. Required back-off.

4. Conclusion

This study has shown that the proposed optimization method can be used for robust multi-objective optimization of chromatographic rare earth element separation, and provided with expected performance changes for the robustification of the studied process. It has been highlighted that the studied system is highly un-robust, and that the system’s lack of robustness is largely due to the neighbouring peaks’ proximity to the product pooling horizon. The system can only cope with slight process disturbances, which in turn demands use of process equipment with high reliability. We show how the optimal solution of a chromatographic separation is affected by introducing robustness in a brute force manner. For future studies, it would be of interest to employ the proposed optimization method on additional chromatography schemes, as well as other chromatography separation applications than rare earth elements.

Acknowledgements

This study has been performed within ProOpt and Process Industry Centre at Lund University, and financed by VINNOVA.

References

[1] Gupta, C.K. and Krishnamurthy, N. (2004) Extractive Metallurgy of Rare Earths. CRC Press, Boca Raton, Florida.

[2] Castor, S.B. and Hedrick, J.B. (2006) Rare Earth Elements. Industrial Minerals Volume, 7th Edition, Society for Mining, Metallurgy, and Exploration, Littleton, Colorado, 769-792.

[3] McGill, I. (2000) Rare Earth Elements. Ullmann’s Encyclopedia of Industrial Chemistry, John Wiley & Sons, Hoboken, New Jersey.

[4] Massari, S. and Ruberti, M. (2013) Rare Earth Elements as Critical Raw Materials: Focus on International Markets and Future Strategies. Resources Policy, 38, 36-43.

[5] Han, K.N., Kellar, J.J., Cross, W.M. and Safarzadeh, S. (2014) Opportunities and Challenges for Treating Rare-Earth Elements. Geosystem Engineering, 17, 178-194.

https://doi.org/10.1080/12269328.2014.958618

[6] Stegen, K.S. (2015) Heavy Rare Earths, Permanent Magnets, and Renewable Energies: An Imminent Crisis. Energy Policy, 79, 1-8.

[7] Golev, A., Scott, M., Erskine, P.D., Ali, S.H. and Ballantyne, G.R. (2014) Rare Earths Supply Chains: Current Status, Constraints and Opportunities. Resources Policy, 41, 52-59.

[8] Nassar, N., Du, X. and Graedel, T. (2015) Criticality of the Rare Earth Elements. Journal of Industrial Ecology, 19, 1044-1054. https://doi.org/10.1111/jiec.12237

[9] Henderson, P. (2013) Rare Earth Element Geochemistry. Elsevier, Amsterdam, Netherlands.

[10] Max-Hansen, M., Knutson, H.K., Jonsson, C., Degerman, M. and Nilsson, B. (2015) Modeling Preparative Chromato-Graphic Separation of Heavy Rare Earth Elements and Optimization of Thulium Purification. Advances in Materials Physics and Chemistry, 5, 151.

https://doi.org/10.4236/ampc.2015.55016

[11] Knutson, H.K., Max-Hansen, M., Jonsson, C., Borg, N. and Nilsson, B. (2014) Experimental Productivity Rate Optimization of Rare Earth Element Separation through Preparative Solid Phase Extraction Chromatography. Journal of Chromatography A, 1348, 47-51.

https://doi.org/10.1016/j.chroma.2014.04.085

[12] Ling, L. and Wang, N.H.L. (2015) Ligand-Assisted Elution Chromatography for Separation of Lanthanides. Journal of Chromatography A, 1389, 28-38.

https://doi.org/10.1016/j.chroma.2015.02.004

[13] Ojala, F., Max-Hansen, M., Kifle, D., Borg, N. and Nilsson, B. (2012) Modelling and Optimisation of Preparative Chromatographic Purification of Europium. Journal of Chromatography A, 1220, 21-25. https://doi.org/10.1016/j.chroma.2011.11.028

[14] Max-Hansen, M., Ojala, F., Kifle, D., Borg, N. and Nilsson, B. (2011) Optimization of preparative Chromatographic Separation of Multiple Rare Earth Elements. Journal of Chromatography A, 1218, 9155-9161.
https://doi.org/10.1016/j.chroma.2011.10.062

[15] Kifle, D., Wibetoe, G., Froseth, M. and Bigelius, J. (2013) Impregnation and Characterization of High Performance Extraction Columns for Separation of Metal Ions. Solvent Extraction and Ion Exchange, 31, 668-682.
https://doi.org/10.1080/07366299.2013.806737

[16] Krattli, M., Muller-Spath, T., Ulmer, N., Strohlein, G. and Morbidelli, M. (2013) Separation of Lanthanides by Continuous Chromatography. Industrial & Engineering Chemistry Research, 52, 8880-8886. https://doi.org/10.1021/ie3031482

[17] Dybczy′ nski, R.S., Kulisa, K., Pyszynska, M. and Bojanowska-Czajka, A. (2015) New Reversed Phase-High Performance Liquid Chromatographic Method for Selective Separation of Yttrium From All Rare Earth Elements Employing Nitrilo-Triacetate Complexes in Anion Exchange Mode. Journal of Chromatography A, 1386, 74-80.
https://doi.org/10.1016/j.chroma.2015.01.091

[18] Knutson, H.K., Holmqvist, A. and Nilsson, B. (2015) Multi-Objective Optimization of Chromatographic Rare Earth Element Separation. Journal of Chromatography A, 1416, 57-63. https://doi.org/10.1016/j.chroma.2015.09.010

[19] Andersson, N., Knutson, H.K., Max-Hansen, M., Borg, N. and Nilsson, B. (2014) Model-Based Comparison of Batch and Continuous Preparative Chromatography in the Separation of Rare Earth Elements. Industrial & Engineering Chemistry Research, 53, 16485-16493.

https://doi.org/10.1021/ie5023223

[20] Jakobsson, N., Degerman, M. and Nilsson, B. (2005) Optimisation and Robustness Analysis of a Hydrophobic Interaction Chromatography Step. Journal of Chromatography A, 1099, 157-166. https://doi.org/10.1016/j.chroma.2005.09.009

[21] Karlsson, D., Jakobsson, N., Axelsson, A. and Nilsson, B. (2004) Model-Based Optimization of a Preparative Ion-Exchange Step for Antibody Purification. Journal of Chromatography A, 1055, 29-39. https://doi.org/10.1016/j.chroma.2004.08.151

[22] Nagrath, D., Bequette, B.W., Cramer, S. and Messac, A. (2005) Multiobjective Optimization Strategies for Linear Gradient Chromatography. AIChE Journal, 51, 511-525.

https://doi.org/10.1002/aic.10459

[23] Gétaz, D., Butté, A. and Morbidelli, M. (2013) Model-Based Design Space Determination of Peptide Chromatographic Purification Processes. Journal of Chromatography A, 1284, 80-87. https://doi.org/10.1016/j.chroma.2013.01.117

[24] Osberghaus, A., Hepbildikler, S., Nath, S., Haindl, M., von Lieres, E. and Hubbuch, J. (2012) Optimizing a Chromatographic Three Component Separation: A Comparison of Mechanistic and Empiric Modeling Approaches. Journal of Chromatography A, 1237, 86-95.

https://doi.org/10.1016/j.chroma.2012.03.029

[25] Sreedhar, B., Wagler, A. and Kaspereit, M., Seidel-Morgenstern, A. (2013) Optimal cut-Times Finding Strategies for Collecting a Target Component from Overloaded Elution Chromatograms. Computers & Chemical Engineering, 49, 158-169.

https://doi.org/10.1016/j.compchemeng.2012.09.009

[26] Ben-Tal, A., El Ghaoui, L. and Nemirovski, A. (2009) Robust Optimization. Princeton University Press, Princeton. https://doi.org/10.1515/9781400831050

[27] Mota, J.P., Araújo, J.M., Rodrigues, R., et al. (2007) Optimal Design of Simulated Moving-Bed Processes under Flow Rate Uncertainty. AIChE Journal, 53, 2630-2642.

https://doi.org/10.1002/aic.11281

[28] Nestola, P., Silva, R.J., Peixoto, C., Alves, P.M., Carrondo, M.J. and Mota, J.P. (2015) Robust Design of Adenovirus Purification by Two-Column, Simulated Moving-Bed, Size-Exclusion Chromatography. Journal of Biotechnology, 213, 109-119.

https://doi.org/10.1016/j.jbiotec.2015.01.030

[29] Holmqvist, A., Andersson, C., Magnusson, F. and Kesson, J. (2015) Methods and Tools for Robust Optimal Control of Batch Chromatographic Separation Processes. Processes, 3, 568-606. https://doi.org/10.3390/pr3030568

[30] Degerman, M., Westerberg, K. and Nilsson, B. (2009) A Model-Based Approach to Determine the Design Space of Preparative Chromatography. Chemical Engineering & Technology, 32, 1195-1202. https://doi.org/10.1002/ceat.200900102

[31] Westerberg, K., Borg, N., Andersson, N. and Nilsson, B. (2012) Supporting Design and Control of a Reversed-Phase Chromatography Step by Mechanistic Modeling. Chemical Engineering & Technology, 35, 169-175.
https://doi.org/10.1002/ceat.201000505

[32] Close, E.J., Salm, J.R., Bracewell, D.G. and Sorensen, E. (2014) A Model Based Approach for Identifying Robust Operating Conditions for Industrial Chromatography with Process Variability. Chemical Engineering Science, 116, 284-295.
https://doi.org/10.1016/j.ces.2014.03.010

[33] Logist, F., Houska, B., Diehl, M. and Van Impe, J.F. (2011) Robust Multi-Objective Optimal Control of Uncertain (Bio) Chemical Processes. Chemical Engineering Science, 66, 4670-4682. https://doi.org/10.1016/j.ces.2014.03.010

[34] Schmidt-Traub, H., Schulte, M. and Seidel-Morgenstern, A. (2012) Preparative Chromatography. 2nd Edition, Wiley-VCH, Weinheim.
https://doi.org/10.1002/9783527649280

[35] Guiochon, G., Felinger, A. and Shirazi, D.G. (2006) Fundamentals of Preparative and Nonlinear Chromatography. Academic Press, New York.

[36] Degerman, M. (2009) Design of Robust Preparative Chromatography. Lund University, Lund.