Multi-objective optimization and decision visualization of batch stirred tank reactor based on spherical catalyst particles

This paper presents a Bayesian approach rooted algorithm oriented to the properties of multi-objective optimization problems. The performance of the developed algorithm is compared with several other multi-objective optimization algorithms. The approach is applied to the multiobjective optimization of a batch stirred tank reactor based on spherical catalyst microreactors. The microbioreactors are computationally modeled by a two-compartment model based on reaction– diffusion equations containing a nonlinear term related to the Michaelis–Menten enzyme kinetics. A two-stage visualization procedure based on the multi-dimensional scaling is proposed and applied for the visualization of trade-off solutions and for the selection of favorable configurations of the bioreactor.


Introduction
Bioreactors, based on microparticles containing immobilized enzyme, permit a specific substrates conversion to a certain product, a use of small volumes of samples and reagents, reduced costs, short processing time and system compactness [18,29].For the development and improvement of highly efficient and productive bioreactors usually a number of physical and biochemical characteristics should be measured and considered [13,18,35].
With the aid of computer tools, the efforts for the design and optimization of bioprocesses and bioelectronic devices can be remarkably reduced [33].The multi-objective optimization of bioprocesses and biosystems has been successfully applied in different applications, particularly, for increasing the productivity and yield of bioreactors [22] and for the optimal design of amperometric biosensors [5].The importance of the multiobjective optimization in chemical and biochemical engineering permanently increases due to the development of new methods sustained by increased computational resources [21,32].
Multi-objective optimization tools provide a mechanism to obtain a certain number of trade-off solutions, known as Pareto-optimal solutions [25].Establishing an efficient approach to find a set of solutions with good trade-off among different objectives has a great practical significance, as these allow engineers to gain insight into the key characteristics of potentially good configurations before moving on to more detailed simulations and pilot commercial tests [16,22,32].Trade-off curves as a visualisation of trade-off solutions are widely used for learning and making decisions when designing products [23].
Multi-objective optimization models are often developed on a basis of mathematical modeling and numerical simulation of the processes [21,22,24,32,33].The selection of multi-objective optimization methods to tackle with those problems depends on the mathematical model of the considered reactor [10,24].For the models of relatively low computational complexity are appropriate metaheuristic algorithms, e.g., genetic algorithms and particle swarm optimization [10,12,17,19,26].
However, the mathematical models of bioreactors based on an array of enzyme-loaded particles are described by transient nonlinear partial differential equations and the numerical solution often requires long lasting computations [2-4, 14, 34].In such cases, the objective functions of the respective optimization problem is so called expensive, and Bayesian approach-based methods can be most appropriate for those optimization problems [24,25].
This paper presents a Bayesian approach rooted algorithm oriented to the properties of multi-objective optimization problems [5].The performance of the developed algorithm is compared with several other multi-objective optimization algorithms.The approach is applied to the multi-objective optimization of a batch stirred tank reactor based on spherical catalyst microreactors [4].The following three objectives were optimized: the reactor operation time and the substrate and enzyme amounts.A two-stage visualization procedure based on the multi-dimensional scaling is proposed and applied for the visualization of trade-off solutions and for the selection of favorable configurations of the bioreactor [37,38].

Mathematical model 2.1 Model domain
A batch stirred tank reactor containing an array of identical spherical microbioreactors placed is considered (see Fig. 1(a)).Assuming a uniform distribution of the microreactors, a unit cell can be considered [3,4].An enzyme-loaded microreactor (MR) occupies the center of the unit cell.The principal structure of the cross section of the modeled unit cell is presented in Fig. 1 In the enzyme-loaded MR, the substrate (S) combines reversibly with an enzyme (E) to form a complex (ES), then the complex dissociates into the product (P) and the enzyme is regenerated [4,35], Since in batch stirred tank reactors the reaction product is usually produced at the same rate as the substrate is consumed, the reactor dynamics can be qualitatively expressed by dynamics only of the substrate concentration [6,13].

Governing equations
Assuming the symmetrical geometry of the MR and homogenized distribution of the enzyme inside the porous MR, the mathematical model can be described in a one-dimensional radial domain [4].Applying the quasi-steady-state approximation of reaction (2), the dynamics of the substrate concentration in MR is described by the following nonlinear reaction-diffusion equation: where s r (r, t) is the concentration of the substrate in the MR, r and t stand for the space and time, respectively, ∆ is the Laplace operator, r 0 is the radius of the MR, D r is the effective diffusion coefficient, e 0 is the total concentration of the enzyme, K M is the Michaelis constant, and k cat is the turnover number [13,35].The product of k cat and e 0 (V max = k cat e 0 ) refer to the maximal enzymatic rate [1,35].
Outside the MR, a thin diffusion shell adjacent to the MR surface remains at a constant thickness h 1 = r 1 − r 0 [35], where s d (r, t) is the concentration of the substrate in the diffusion shell, and D d is the corresponding diffusion coefficient.
The substrate is assumed to be uniformly distributed throughout the outside of the diffusion shell and its concentration depends only on time [7], , t > 0, where s b (t) is the substrate concentrations in the convective shell, q is the ratio of the volume of the convective enclosure (r 1 r r 2 ) to the area of the outer surface of the diffusion shell (r = r 1 ), The thickness h 2 = r 2 − r 1 of the convective shell indirectly corresponds to the density as well as to the number of the microreactors placed in the container.

Initial and boundary conditions
The process starts (t = 0) when substrate of concentration s 0 is injected into the reactor, s r (r, 0) = 0, 0 r r 0 , The zero-flux boundary condition is used for the MR centre (t > 0), The formal partition coefficient φ is used in the matching conditions to describe the specificity in the concentration distribution, and the flux of the substrate through the diffusion shell is assumed to be equal to the flux entering the MR surface (t > 0) [15,34], The continuity of the substrate concentration is defined on the boundary between the diffusion and convective shells (r = r 1 , t > 0), http://www.journals.vu.lt/nonlinear-analysis 3 Multi-objective optimization problem

Problem statement
A minimization of time-cost is often sought by designers of biotechnological processes [18,31].The batch time required to achieve a certain conversion of the reactants is usually assumed as the main characteristic of the process duration [13].The reactor should convert as much substrate as possible within the shortest possible time.
In some applications of bioreactors, enzymes are available only in microgram to milligram quantities and are very expensive [13,35].In such applications, the minimization of the enzyme usage is of crucial importance.In the case of biocatalytic microreactor shown in Fig. 1(b), the amount of the enzyme loaded into the MR is equal to the product of the initial enzyme concentration e 0 and the microreactor volume (V r = 4πr 3 0 /3).The amount of enzyme loaded into all the microreactors placed into the container equals e 0 V r N , where N is the number of microreactors.
For the optimization, without losing generality, the volumetric density N V of microreators placed in the container (N V = N/V ) can be used instead of the total number N = 3V /(4πr 3  2 ) of microreactors.Thus, the amount of enzyme used in N V microreactors (used per volume unit of the container) equals 2 ) = e 0 r 3 0 /r 3 2 .Accordingly, the amount of the substrate in the bulk per volume unit of the container can be expressed as follows: 1 )/r 3 2 .The optimal design of the batch reactor mathematically is stated as a multi-objective optimization problem with three objective functions, ϕ 1 (r 0 , e 0 , s 0 , r 2 ) = t 0.9 = t: s b (t) = 0.1s 0 , where ϕ 1 (•) stands for the batch time t 0.9 required to convert 90% of the initial amount of the substrate [4], ϕ 2 (•) is the amount of the substrate per volume unit converted to the product, and ϕ 3 (•) is the total amount of the enzyme used per volume unit of the reactor.The first and third objectives should be minimized while the second one should be maximized.The appropriate intervals of the decision variables for the optimal design problem are given in Table 1.Assuming the highly stirred reactor, the thickness h 1 of the diffusion shell can be assumed as a constant parameter [3,4].The radius r 2 of the unit cell can be expressed via (independent) decision variables r 0 , h 2 , and the parameter h 1 as follows: Additionally, only reactor configurations with reasonable batch time were considered.If the time required to convert 90% of the initial amount of the substrate is less than the time limit T max = 10 4 seconds (i.e., t 0.9 10 4 s), the reactor configuration was considered as acceptable [4,13,35].Otherwise the problem solution was excluded from a further analysis.
The following typical values of the model parameters were kept constant [1,11,13,35]:

Optimization algorithm
A main challenge of the three-objective ( 13) optimization problem is the computational complexity of the objective function ϕ 1 (x), the computation of which includes the numerical integration of the transient nonlinear governing equations [4].Because of the nonlinearity of the initial boundary value problem ( 3)-( 12), the bioreactor action was simulated numerically using the finite difference method with explicit scheme [3,4].
Although explicit difference schemes have the strict stability limitations, these schemes have a convenient algorithm of the calculation and are simple for programming [8].The computation of a single value of ϕ 1 (x) using computer with Intel Xeon X5650 2.66GHz processor can take up to 10 minutes.
For the solution of applied black-box multi-objective optimization prevail meta-heuristic methods [12,30].However, the objective functions of the considered problem are computationally expensive, and the solution of this problem by a meta-heuristic method would take prohibitive long computing time.For such, so called expensive, objective functions, Bayesian methods proved to be efficient [25].On the other hand, the inherent complexity of standard implementations of Bayesian methods limits the number of values of the objective functions which can be processed.The modest number of iterations seems sufficient for the solution of various applied problems, e.g., up to 60 evaluations of values of the objective functions in [24], 100 function evaluations in [27], and 200 functions evaluations in [20].However, such a number of function values is not sufficient for the appropriate representation of the Pareto front of the considered three-objective optimization problem.Therefore, a partition based implementation of the P-algorithm was implemented.The P-algorithm is a Bayesian approach rooted algorithm defined by the maximum probable improvements at each iteration.The improvement means a correction of the approximation of the Pareto front, and the improvement probability is computed with respect to the accepted statistical model of objective functions.The developed implementation is based on the theoretical multi-objective method proposed in [36] and adapts the implementation of the single-objective P-algorithm proposed in [9].The implementation of the P-algorithm in [9] is based on the rectangular partition of the feasible region, http://www.journals.vu.lt/nonlinear-analysiswhere a hyper-rectangle for the bisection is selected according to maximum probability of the improvement of the estimate of global minimum.In the newly developed multiobjective generalization, the selection criterion is the probability of the correction of the approximation of the Pareto front.The computational burden of the proposed algorithm is essentially lower than that of standard implementations of Bayesian algorithms.The further reduction of computational burden is planned by the replacing the rectangular partition with the simplicial partition; see, e.g., [39].
The optimization algorithm, as well as the bioreactor simulation algorithm, was implemented in C++ [4,5,28].The computing time referenced below is measured experimenting using computer with Intel Xeon X5650 2.66GHz processor.

Numerical example
The following bi-objective minimization problem is solved to illustrate the performance of the developed algorithm by comparison its testing results with the testing results of several other multi-objective optimization algorithms: The problem is considered hard because of non-convexity, and of the coincidence of the set of Pareto optimal decisions with the boundary of the feasible region.Let us note that the Pareto optimal solution (1, 0) T corresponds to the set of optimal decisions P D = {x: 0 x   The performance metrics of the proposed algorithm are presented in Table 2, where the performance of metrics of several other algorithms are also presented for comparison [25].The following notation is used in Table 2: NFE-number of evaluations of the objective functions, ε -tolerance of the termination condition, NGEN-number of generation in the termination criterion of the genetic algorithm, HV-hypervolume computed with respect to the reference point (1, 1) T , NP-number of found non-dominated solutions, UDuniformity of distribution of points approximating the Pareto front.For the convenience of readers the definition of the metric UD follows [25]: where k is the number of found non-dominated points, and d i is the minimum distance from the point of the set of non-dominated points indexed by i to the other point of this set, i = 1, . . ., k.

Solution of the three-objective optimization problem
The considered optimization problem is three-objective (13) with four decision variables (Table 1).The variables in the problem description were defined in physical units.The equations of the mathematical model are well conceivable when variables are presented in natural units, however, the feasible values of variable in this case differ in several orders.Since the transition to the logarithmic scale facilitates the proper partitioning of the feasible region by the considered optimization algorithm, the variables should be rescaled.The new variables are: x = (x 1 , . . ., x 4 ) T , x 1 = lg r 0 , x 2 = lg e 0 , x 3 = lg s 0 , and x 4 = lg h 2 , and the potential feasible region is The optimization algorithm uses internal scales, where the variables vary in the interval [0, 1].We used the term "potential" in the definition of A since some infeasible subsets of A remain not defined explicitly.The infeasibility is determined during the simulation of the bioreactor action, meaning that the reactions in bioreactor are not completed during the time limit T max , i.e., t 0.9 > T max .Since a returned undefined function value can crash optimization process, the algorithm computing the objective functions was stopped by reaching the predefined bioreactor batch time t 0.9 = T max and returned the values of the objective functions corresponding to the maximum feasible simulation time.
http://www.journals.vu.lt/nonlinear-analysisWe expected to obtain reasonable approximation of the Pareto front after 1000 evaluations of the objective functions.The computations continued 4342 minutes.However, 351 computations out of thousand were indicated as infeasible.The number of non-dominated points was equal to 124.The corresponding representation of the Pareto front is shown in Fig. 3.
The general shape of Pareto front shown in Fig. 3 is not very helpful for making a decision.An interactive extension of the optimization algorithm also does not seem promising because of long lasting computations of the objective functions.
The considered optimization problem is difficult partly due to the implicitly infeasible points.The description of a region of the infeasible points could increase the optimization efficiency by avoiding the redundant time consuming computations.Methods of the visualization of multi-dimensional data were shown helpful in the solution of similar problems [5,37,38].
The location of the infeasibility region in the hypercube can be understood by the analysis of the two-dimensional images of the infeasible points with respect to the image of the vertices of the hypercube.Therefore, a visualization method is needed to produce an image with clearly expressed structural properties of the hypercube.The multi-dimensional scaling (MDS) based visualization meets these requirements [37,38].However, the vertices constitute only a small fraction of the visualization data (351 infeasible points, 124 Pareto optimal solutions, and 16 vertices), therefore, the original structure of the set of vertices is not hold in the image of whole set of data.Since we cared about the position of the points in question with respect to vertices but not about the mutual distances between them, we have applied a two-stage visualization procedure.First, an image of the set of vertices complemented with the average of the infeasible vectors was computed by the hybrid algorithm of MDS [38].The sequential MATLAB implementation of the MDS algorithm was used as suitable for the modest dimensionality of the data.In the second stage, the image z of an infeasible point (or an optimal Pareto Let us note the closeness of the images of the infeasible points (red crosses in Fig. 4) to the vertices (indicated by the left pointing triangles) the third coordinate of which is equal to −8.Based on this observation we can guess that the infeasibility is related to the small enzyme concentration.Contrary, the images of the Pareto optimal solutions (blue plus signs) are closer to the images of vertices, the third coordinate of which is equal to −4.Consequently, we can guess that the enzyme concentration of the majority of Pareto optimal decisions is distant from the lower bound of the enzyme concentration given in Table 1, and that only a few Pareto optimal decisions are closely located to the infeasible points.

Reformulation of the problem
Graphical presentation of the Pareto front of bi-objective optimization problems excellently aids decision makers to compare quantitatively the available alternative before http://www.journals.vu.lt/nonlinear-analysis decision making.However, in the case of three objective problems, similar 3D graphs are useful rather for a qualitative illustration of the Pareto front than for the tradeoff between potential decisions.Fortunately, the considered problem ( 13) and ( 16) can be reduced to several bi-objective problems.The structure of the considered optimization problem is favorable to reformulation which not only reduces the number of objectives but also the number of variables.The idea is following.Let fix several favorable values of third objective function, and compute and draw the Pareto fronts of the modified, bi-objective, optimization problems.An appropriate solution can be found by interpolating between the drawn fronts.
Let us reformulate the original optimization problem by introducing a new independent variable r 2 .h 2 becomes dependent and is expressed through variables r 2 , r 0 and constant parameter h 1 , h 2 = r 2 − r 0 − h 1 .Instead of functions φ(•), we will consider their (decimal) logarithms where x 1 = lg r 0 , x 2 = lg e 0 , x 3 = lg s 0 and x 4 = lg r 2 .
The feasible region is expressed by the bounds defined for the original independent variables, Let us fix an appropriate value of f 3 (•), E = f 3 (x).The last equation of ( 17) can be rewritten in the form meaning that the number of independent variables can be reduced to three, and the optimization problem (17) can be reduced to the following bi-objective parametric minimization problem: where the feasible region X of three independent variables is defined as follows:

06) .
Nonlinear Anal.Model.Control, 24(6):1019-1033 The reformulated problem is more convenient for an analysis since the number of variables and objectives is smaller than in the original problem.The time-consuming computation of the first objective f 1 (•) remains the main difficulty of the problem.However, an acceptable solution can be achieved with smaller number of the evaluations of the objective functions because of the reduced dimensionality.

Analysis of optimal solutions
Three multi-objective optimization problem (18) was solved with the following three values of the parameter E: −7, −6.5, −6, where E denotes a fixed value of f 3 (•).These values of E corresponds to the following relatively high values of the enzyme amount ϕ 3 required per volume unit of the reactor: 10 −7 , 10 −6.5 and 10 −6 M.
The multi-objective optimization algorithm used for the three objective optimization problem (13) was also used for the reduced optimization problem (18).The termination condition was defined by the budget of computations of the objective functions equal to 500.The obtained Pareto fronts are presented in Fig. 5.One can see in Fig. 5 not fully smooth curves since the global optimization algorithm stops not precisely reaching the Pareto front.The specifications would be possible using a local minimization algorithm, however, this would require long lasting computations needed for the simulation of the reactor action.
The optimization results presented in Fig. 5 can be used for the evaluation of a potential configuration of a bioreactor (trade-off solution).For example, a bioreactor with the following parameters could be of interest: the batch time ϕ 1 = t 0.9 = 3600 s, the amount of the reduced substrate ϕ The other parameters of these two solutions are presented in the first two lines of Table 3.Since the desirable value of ϕ 2 (x) is close to the average of the corresponding values at the considered Pareto fronts it can be expected that the average of input parameters of these two versions of bioreactors appropriately approximates the desirable version.
The average values of the parameters are presented in the third line of Table 3, and the computed values of the objective functions are presented in the fourth line of this table.Indeed, the computed values of the objectives are quite close to the desirable ones.The neighbourhood of the obtained solution can be further analyzed either by the methods of visualization or specified by a local optimization method.The Pareto fronts computed for the reformulated two objective optimization problem with several fixed values of f 3 (•) facilitate the analysis of trade-off solutions and aids a designer to determine an optimal configuration of a bioreactor.

Conclusions
A multi-objective optimization can be successfully used for the development of efficient and productive batch stirred tank reactors based on catalyst particles.The development and improvement of such bioreactors involve consideration of the simultaneous optimization of several objectives, some of whose are conflicting.The mathematical and corresponding numerical models of a bioreactor can be used for the evaluation of the objectives.
The stated problem of the multi-objective optimization of bioreactors is difficult to solve due to relatively high computational complexity of the objectives as the solutions of transient nonlinear reaction-diffusion problems.A Bayesian approach rooted algorithm defined by the maximum probable improvements at each iteration is most appropriate for solving the stated optimization problem (Table 2).A two-stage visualization procedure based on the multi-dimensional scaling is a suitable method for the visualization of tradeoff solutions and the selection of favorable configurations of the bioreactor (Fig. 4).
For clearer analysis and graphical representation of the trade-off solutions, the multiobjective optimization problem can be reduced to several bi-objective problems with several fixed values of a certain objective, and the interpolation can be used for calculating an appropriate solution (Fig. 5, Table 3).

Figure 1 .
Figure 1.Principal structure of a container containing spherical microreactors (a) and the cross section of a modeled unit cell (b).

1 1 ,
x 2 = 0}.The computed approximation of the Pareto front and the corresponding set of decisions are shown in Fig. 2. In the right-hand side of Fig. 2, the rectangular partition of the feasible region is shown at the vertices of which are computed the values of the objective functions.

Figure 2 .
Figure 2. Spaces of objectives (a) and variables (b) illustrating the results optimizing the problem (15).

Figure 3 .
Figure 3.The representation of Pareto front.

Figure 4 . 16 i=1t
Figure 4. Visualization of the distribution of infeasible points (indicated by •), Pareto optimal solutions (+) and vertices of the hypercube ( ).solution) x was defined as follows:

Table 2 .
Summary of the results of numerical examples.

Table 3 .
Data and results of a candidate trade-off solution.