ABSTRACT

Rare earth elements are strategic commodities in many countries, and an important resource for the growing modern technology industry. As such, there is an increasing interest for development of rare earth element processing, and this work is a part of further development of chromatography as a rare earth element separation process method. Process optimization is pivotal for process development, and it is common that several competing objectives must be regarded. Chromatographic separation processes often consider competing objectives, such as productivity, yield, pool concentration and modifier consumption, which leads to Pareto optimal solutions. Adding robustness to a process is of great importance to account for process disturbances and uncertainties but generally comes with reduced performance of the other process objectives as a trade off. In this study, a model-based robust multi-objective optimization was carried out for batch-wise chromatographic separation of the rare earth elements samarium, europium and gadolinium, which was considered highly un-robust due to the neighbouring peaks proximity to the product pooling horizon. The results from the robust optimization were used to chart the required operation point changes for keeping the amount of failed batches at an acceptable level when a certain level of process disturbance was introduced. The loss of process performance due to the gained robustness was found to be in the range of 10% - 20% reduced productivity when comparing the robust and un-robust Pareto solutions at Pareto points with identical yield. The methodology presented shows how to increase robustness to a highly un-robust system while still keeping multiple objectives at their optima.

Keywords:

Rare Earth Elements, Chromatography, Multi-Objective Optimization,

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}_{\max ,\alpha }\left[1-\sum _{\gamma \in \left\{\text{Sm,Eu,Gd}\right\}}\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,α}}$ 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,α}}\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{if}\alpha \in \left\{\text{Sm,Eu,Gd}\right\}\mathrm{,}\hfill \\ c_{\text{mix}\mathrm{,}S}{v}_{\text{int}}\hfill & \text{if}\alpha =S\mathrm{,}\hfill \end{array}\end{array}$ (3)

$\frac{\partial c_{\alpha }}{\partial z}\left(t\mathrm{,}{z}_f\right)=\mathrm{0,}\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{if}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{if}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

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=2V_c{\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)}{\sum _{\beta \in \left\{\text{Sm},\text{Eu},\text{Gd}\right\}}{\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{<mo\; stretchy="false">\; (\; </mo>}t\mathrm{,}{z}_k\mathrm{<mo\; stretchy="false">\; )\; </mo>,}{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{.}-\left(\omega \underset{t_0}{\overset{t_f}{\int }}\frac{\text{d}{P}_{\text{Eu}}}{\text{d}t}\text{d}t+\left(1-\omega \right)\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{.}p=\left(\Delta t_{\text{load}}\mathrm{,}{u}_0\mathrm{,}{u}_f\right)\in {ℝ}^3\mathrm{,}$

$\text{s}\text{.t}\text{.}{p}_L\le p\le p_U\mathrm{,}$ (10b)

$\left(x\mathrm{,}{t}_c\mathrm{,}{t}_f\right)=\text{argmin}\text{.}-\left(\omega \underset{t_0}{\overset{t_f}{\int }}\frac{\text{d}{P}_{\text{Eu}}}{\text{d}t}\text{d}t+\left(1-\omega \right)\underset{t_0}{\overset{t_f}{\int }}\frac{\text{d}{Y}_{\text{Eu}}}{\text{d}t}\text{d}t\right)\mathrm{,}$ (10c)

$\text{w}\text{.r}\text{.t}\text{.}\left(t_c\mathrm{,}{t}_f\right)\in {ℝ}^2\mathrm{,}$

$\text{s}\text{.t}\text{.}\dot{x}=F\left(t\mathrm{,}x\left(t\right)\mathrm{,}{t}_c\mathrm{,}{t}_f\mathrm{,}p\right)\mathrm{,}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{,}{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{,}\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, $\tilde{p}$ , on the free operating parameters, $p$ (i.e $\Delta t_{\text{load}}\mathrm{,}{u}_0$ and $u_f$ ), and define ${\tilde{X}}_{\text{Eu}}$ as the cumulative purity distribution of the model responses that are produced from $\tilde{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{.}\underset{-\infty }{\overset{X_{\text{Eu}\mathrm{,}L}+X_{\text{BF}}}{\int }}{\tilde{X}}_{\text{Eu}}\text{d}{X}_{\text{Eu}}-{\Phi }_{X_{\text{Eu}}}\mathrm{,}$ (11a)

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

$\text{s}\text{.t}\text{.}\dot{\tilde{x}}=\tilde{F}\left(t\mathrm{,}x\left(t\right)\mathrm{,}{t}_c\mathrm{,}{t}_f\mathrm{,}\tilde{p}\right)\mathrm{,}$ (11b)

$\tilde{p}~N\left(p\mathrm{,}{\sigma }_p^2\right)\mathrm{,}$ (11c)

$p=\text{argmin}\text{.}-\left(\omega \underset{t_0}{\overset{t_f}{\int }}\frac{\text{d}{P}_{\text{Eu}}}{\text{d}t}\text{d}t+\left(1-\omega \right)\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{.}p=\left(\Delta t_{\text{load}}\mathrm{,}{u}_0\mathrm{,}{u}_f\right)\in {ℝ}^3\mathrm{,}$

$\text{s}\text{.t}\text{.}{p}_L\le p\le p_U\mathrm{,}$ (11e)

$\left(x\mathrm{,}{t}_c\mathrm{,}{t}_f\right)=\text{argmin}\text{.}-\left(\omega \underset{t_0}{\overset{t_f}{\int }}\frac{\text{d}{P}_{\text{Eu}}}{\text{d}t}\text{d}t+\left(1-\omega \right)\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{.}\left(t_c\mathrm{,}{t}_f\right)\in {ℝ}^2\mathrm{,}$

$\text{s}\text{.t}\text{.}\dot{x}=F\left(t\mathrm{,}x\left(t\right)\mathrm{,}{t}_c\mathrm{,}{t}_f\mathrm{,}p\right)\mathrm{,}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{,}{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{,}\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, $\tilde{x}$ , for an uncertainty set of the free operating parameters, $\tilde{p}$ , to evaluate the cumulative distribution function of ${\tilde{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 ${\tilde{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

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, $\tilde{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, ${\tilde{X}}_{\text{Eu}}$ , of the uncertainty set.

Then, an initial investigation of the back-off terms impact on ${\tilde{X}}_{\text{Eu}}$ was conducted by creating new Pareto fronts with an incrementally increased back-off and observing how ${\tilde{X}}_{\text{Eu}}$ changes when $\tilde{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 ${\tilde{X}}_{\text{Eu}}$ based on $\tilde{x}$ and $\tilde{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 $\tilde{x}$ , $\tilde{p}$ and cumulative distribution functions of ${\tilde{X}}_{\text{Eu}}$ to be evaluated. This continues until a $X_{\text{BF}}$ that produces a cumulative distribution function of ${\tilde{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.

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.

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

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.

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.

Cite this paper

Knutson, H.-K., Holmqvist, A., Andersson, N. and Nilsson, B. (2017) Robust Multi-Objective Optimization of Chromatographic Rare Earth Element Separation. Advances in Chemical Engineering and Science, 7, 477-493. https://doi.org/10.4236/aces.2017.74034

References