Robust Multi-Objective Optimization for Response Surface Models Applied to Direct Low-Value Natural Gas Conversion Processes

The high proportion of CO2/CH4 in low aggregated value natural gas compositions can be used strategically and intelligently to produce more hydrocarbons through oxidative methane coupling (OCM). The main goal of this study was to optimize direct low-value natural gas conversion via CO2-OCM on metal oxide catalysts using robust multi-objective optimization based on an entropic measure to choose the most preferred Pareto optimal point as the problem’s final solution. The responses of CH4 conversion, C2 selectivity, and C2 yield are modeled using the response surface methodology. In this methodology, decision variables, e.g., the CO2/CH4 ratio, reactor temperature, wt.% CaO and wt.% MnO in ceria catalyst, are all employed. The Pareto optimal solution was obtained via the following combination of process parameters: CO2/CH4 ratio = 2.50, reactor temperature = 1179.5 K, wt.% CaO in ceria catalyst = 17.2%, wt.% MnO in ceria catalyst = 6.0%. By using the optimal weighting strategy w1 = 0.2602, w2 = 0.3203, w3 = 0.4295, the simultaneous optimal values for the objective functions were: CH4 conversion = 8.806%, C2 selectivity = 51.468%, C2 yield = 3.275%. Finally, an entropic measure used as a decision-making criterion was found to be useful in mapping the regions of minimal variation among the Pareto optimal responses and the results obtained, and this demonstrates that the optimization weights exert influence on the forecast variation of the obtained response.


Introduction
A mixture of various hydrocarbons such as methane (CH 4 ), ethane (C 2 H 6 ), propane (C 3 H 8 ), butane (C 4 H 10 ) and inert diluents such as molecular nitrogen (N 2 ) and carbon dioxide (CO 2 ) are found in natural gas compositions. Natural gas composition variations are affected by several parameters such as the geographical source, the time of year, and the treatment applied during production or transportation [1][2][3].
CH 4 and CO 2 are believed to contribute to the greenhouse effect [4]. Although the amount of CH 4 in the atmosphere is less than CO 2 , its global warming potential is approximately 25 times greater [5]. Large amounts of methane are widely available in nature in the form of natural gas, although methane is a greatly underutilized resource for chemical and liquid fuels [6]. Large amounts of natural gas reserves are located in remote areas [7,8] and pipelines may not be available to transport that remote gas to regions where it can be used, and liquefication for transportation via ship is quite expensive. Almost 11% of this gas is reinjected, while another 4% is flared or vented [5,9,10], generating waste [10].
The conversion of CH 4 and CO 2 into value-added chemicals has attracted attention from the academic community and industries [4,[11][12][13]. The abundance of these two gases makes them stand out as raw materials for fuels and chemical synthesis. Furthermore, CH 4 is considered to be the cheapest carbon source for the petrochemical industry [14].
Among these reserves, there are large amounts of low-value natural gas containing a high concentration of CO 2 . Currently this natural gas is sold as liquefied natural gas after a carbon dioxide clean up using CO 2 separation facilities. The impurities from both the separation stages are then injected underground to dispose of the waste using a sub-surface aquifer, but at an estimated 49% of the total cost of the project [15]. The high CO 2 /CH 4 ratio in low-value natural gas compositions could be strategically utilized to produce value-added chemicals, such as higher hydrocarbons and liquid fuels without having to separate the CO 2 first [16]. The synthesis of liquid fuels and commodity chemicals from CO 2 is a promising approach for clean energy production. For this reason, much academic and industrial effort has been devoted to exploring efficient means of reducing CO 2 [17].
Several studies have been undertaken with the objective of making conversion of methane viable, either by direct or indirect routes. The indirect routes focus on steam reform to produce synthesis gas (CO + H 2 ), which can be converted into the liquid fuels [18]. The direct route, on the other hand, converts methane into higher hydrocarbons in one step through oxidative methane coupling reactions (OCM). Therefore, it is considered more economically viable and, consequently, it has been the subject of several studies [4,12,14,16,[19][20][21][22][23][24][25][26][27][28][29][30][31]. The general reactions for the formation of C 2 hydrocarbons from methane and carbon dioxide are expressed by [24]: Although the OCM reaction is a highly exothermic reaction, it requires high temperatures and a suitable catalyst because the binding energy of hydrogen-carbon in methane is very large [32] and the bond dissociation energy of CO-O is also high. Because of the low reactivity of CO 2 , the product (C 2 H 4 , C 2 H 6 ) of the OCM reaction, with CO 2 as an oxidant, is less likely to react with CO 2 . Therefore, high C 2 (C 2 H 4 , C 2 H 6 ) selectivity is expected in the CO 2 -OCM [31].
Nevertheless, a possible increase in temperature will result in total oxidation instead of partial oxidation of C 2 hydrocarbons, such as ethylene [6]. We therefore see the importance of catalyst selection for the OCM reaction. The catalyst determines at what temperature and composition the maximum yield and selectivity for the reaction will be defined [30].
Of the various catalysts for CO 2 activation, ceria (CeO 2 ) is attracting increased interest due to its high oxygen storage capacity [11]. Indeed, oxidation and reduction reactions of Ce 4+ and Ce 3+ are effective in activating carbon dioxide to form oxygen active species, while the C 2 selectivity is related to the basicity of a catalyst due to enhanced CO 2 chemisorption on the catalyst surface [11,27]. Wang et al. [25] and Wang et al. [26] proposed using a ceria catalyst modified with CaO. According to these authors the CaO/CeO 2 is the most effective catalyst for conversion at high temperatures of CH 4 to C 2 H 6 and C 2 H 4 by CO 2 among a series of CeO 2 catalysts modified with alkali and alkaline earth metal oxides. Istadi and Amin [28] and Istadi and Amin [16] proposed using CaO/CeO 2 and MnO/CeO 2 catalysts in the multi-objective optimization of CO 2 -OCM process. According to these authors, the optimal operating parameters, such as the CO 2 /CH 4 ratio and reactor temperature, and the catalyst compositions in the ceria catalyst, provide essential information for industrial CO 2 -OCM processes.
At the same time, multi-objective optimization allows one to find the point that represents the final solution to the problem, in addition to treating the reliability of that solution as an important factor. It then follows that the variance of the prediction would be of great concern [33].
The aim of this study is to optimize the conversion of natural gas, rich in carbon dioxide, using robust multi-criteria decision-making processes, based on an entropic measure, to determine the ideal Pareto point as the final solution to the problem. The normal boundary intersection method (NBI) together with the mixture design of experiments (MDE), will be used to optimize these responses simultaneously.

Multi-Objective Optimization
Industrial processes, when translated into optimization problems, are usually treated as multi-objective problems, since they involve more than one desirable resource. When the objective functions are not in conflict, each objective function reaches its optimal value, and a solution is found, without the need for any special method [34]. These objectives are often conflicting [35]. A multi-objective optimization (MOP) problem can be formulated in order to study the tradeoffs between these conflicting objectives, as in: where Λ is the vector of the objective functions that comprise the conflicting k criteria and F i is the vector of the decision variables x belonging to the viable set Ω.
Their restrictions are represented in the form of inequalities or equalities, according to: where g r and h q are the inequality and equality constraint functions, respectively and I and J are the index sets containing as many elements as there are inequality and equality constraints, respectively. When looking for an MOP solution, the goal is to find efficient solutions. Some variations in the concept of efficiency are concepts of local efficiency, weak efficiency, and weakly local efficiency. Hence, a solution x * ∈ X is locally efficient if δ > 0 such that However, for a given decision vector x * to be the solution to an optimization problem, it must satisfy some conditions, which we call optimality conditions. The simplest situation for which we can define the optimality conditions is that in which we have a single function for which we wish to find the optimal point, and there are no constraints. In the case of unconstrained optimization, the necessary first-order condition is: If x * is to be a local minimizer of the function f (x), differentiable at x * , then ∇ f (x * ) = 0. At this point, Miettinen [34] is used, which shows that a function f i : R n → R is differentiable at x * if: . . , n), i.e., all the gradient components, are continuous at x * . Still, with respect to the unconstrained optimization, the necessary second-order condition is: If x * is to be a local minimizer of the function f (x), twice differentiable at x * , then ∇ f (x * ) = 0 and ∇ 2 f (x * ) must be positive semidefinite, i.e., it has eigenvalues (λ i ) greater than or equal to zero. According to Miettinen [34], a function f i : R n → R is twice differentiable at x * if: where: ∇ f i (x * ) is the gradient, the symmetric matrix, n × n; The Hessian matrix of a twice differentiable function consists of second-order partial . . , n), and it can be presented as: .
Furthermore, f i is twice continuously differentiable at x * if all of its second-order partial derivatives are continuous at x * .
With respect to optimality conditions in unconstrained optimization, we see that only the presented conditions are necessary, since the first and second order terms can be null, still leaving doubt about the nature of x * . Therefore, a sufficient condition for x * to be a strict local minimizer of the function f (x), twice differentiable at x * , is ∇ f (x * ) = 0 and ∇ 2 f (x * ) as a positive definite, i.e., it has eigenvalues (λ i ) greater than zero.
To analyze the critical points the second order derivatives of the function must exist and be different from zero, which can be verified by performing a Taylor series expansion around the optimal point. Since the first derivative is null, we have: The values of the function near the critical point depend on the Hessian. In summary, we have the following relation between the Hessian matrix eigenvalue signals and the critical point: This analysis still allows us to deduce the convexity of the function. Similarly, when working with Response Surface Methodology (RSM), the determination of the convexity of a function is done by characterizing the nature of the stationary point. The stationary point is the level of x 1 , x 2 , . . . , x k , that optimizes the predicted response. This point, if it exists, will be the set of x 1 , x 2 , . . . , x k , for which the partial derivatives are equal to zero. A general mathematical solution for locating the stationary point may be obtained. The second-order model is expressed in the matrix notation as [36]: where: The derivative ofŷ with respect to the elements of the vector x equated to zero is: The stationary point described in Equation (12) is the solution of Equation (11): The predicted response at the stationary point is given by: In general terms, the nature of a stationary point can be determined from the sign of the eigenvalues or root characteristics of a given matrix B. The eigenvalues (λ i ) of the matrix B are the solutions to the following equation: Then, these are defined as follows: • if the values of λ i are all negative, the function is concave and x S is a maximum point; • if the values of λ i are all positive, the function is convex, and x S is a minimum point; • if the values of λ i present different signs, the function is neither concave nor convex, and x S is a saddle point.
In multi-objective optimization two approaches can be used to aid in solving problems: (i) converting all objective functions into a single problem; (ii) optimizing one objective, considering other objectives as constraints [37]. In this second approach, the objective function is prioritized and the relevance of the others is relegated to a lower extent.
According to Shahraki and Noorossana [37], the methods that result in obtaining a set of optimal Pareto solutions are recommended, since they provide the best solutions among all other given options in terms of efficiency. The weighted sum method (WSM) is widely used to obtain optimal Pareto solutions in MOP, as it is easy to implement and interpret. However, if the set of Pareto optimal solutions is nonconvex, the frontier becomes nonconvex and discontinuous, forming clusters of solutions in regions of great curvature, yet discontinuous in the solution space. In such situations, the WSM, which is the standard method for generating the Pareto set in MOP, barely finds solutions in the nonconvex section. Moreover, the WSM cannot generate an equally spaced frontier, even if the distribution of weights is uniform [38,39], which can confuse the decision maker by not clarifying the conflicting behavior and the trade-offs between different objective functions.
Das and Dennis [40] presented the Normal Border Intersection (NBI) method as an option that can overcome the disadvantages presented by the WSM method. This method presents the Pareto surface distributed evenly independent of relative scales and the convexity of the objective functions.
However, Das and Dennis [40] argue that a disadvantage, inherent to methods that seek to find a large number of efficient points in MOP, is that these methods cannot find globally Pareto optimal points. The points generated by the NBI are only locally guaranteed as Pareto optimal points. Nevertheless, we will use the NBI method in this study, given its robustness when working with nonconvex problems.
Thus, the NBI method is used to solve the MOP using the following equation [40]: where w is the convex weighting; D is the distance between the Utopia line and the Pareto frontier; F(x) is the vector containing the individual values of the normalized objectives in each run; e is a column vector of one and Φ and Φ are the payoff and normalized payoff matrices, respectively, and can be written as:

Criteria for Defining the Ideal Pareto Solution
When solving an MOP, there are usually an infinite number of efficient solutions that form the ideal set of Pareto (called efficient set) points [34]. The process that seeks to generate optimal Pareto alternatives is known as multi-objective optimization. From a mathematical standpoint, every Pareto ideal point is an acceptable solution for a MOP, if the objective is to obtain a point as the final solution [41].
It is difficult to define the degree of importance to be attributed to each objective [42]. The definition of the weights for each function can be influenced by the preferences of the decision maker. This affects the influence of weights used to determine the relative importance of the functions, in order to identify the most important parameters during the optimization process, and the preferences are selected [43].
The priority given to the criteria significantly affects the final result, since this result depends on the importance attached to each objective [44][45][46]. This can become a problem as decision makers are often unsure about the exact weights of objective functions [44]. Considerable cognitive effort is necessary [47] in order to obtain information on direct preference from the analyst, The Pareto set includes rational options, from which the analyst can select the final solution by comparing several objectives with each other [44]. As such, the search is for a set of ideal solutions in the broadest sense (Pareto optimal). Several techniques from literature address optimal Pareto solutions in the solution space. However, the disadvantage of these methods is the variety of solutions from which one must choose. We therefore can identify a need to bridge the gap between exclusive solutions and Pareto ideal sets [44].
When solving a linear multi-objective optimization problem, Zeleny [48,49] sought to answer the following questions: (i) what is the most preferred solution among the generated, non-dominated and extreme solutions? (ii) Can the set of non-dominated solutions be reduced to consist of fewer points to determine a final decision? To answer these questions, he proposed the "traditional measure of entropy" as a parameter to assess the importance of functions and define the weights to be used in solving the problem.
Following a different approach, some authors (see [50][51][52]) used the Shannon entropy index [53] associated with an error measure, to determine the most preferred Pareto ideal point in a MOP in a vertical turning process, resolved using the NBI method. The authors state that Shannon's entropy index can provide a more reliable assessment of the relative weights of objectives in the absence of analyst preferences. Furthermore, the authors state that when combined with an error measure, it minimized the error of Pareto's preferred point related to individual optimal responses. The weighting metric ξ is obtained by [50,52]: where w i are the weights to be assigned to the objectives that are to be optimized. The Entropy equation can be calculated by [53]: The GPE equation is represented by [54]: where y * i is the value of the Pareto-optimal responses; T i is equal to the defined target and m is equal to the number of objectives.
As discussed earlier, many of the weighting strategies used during the optimization and decision-making process rely on inaccurate and subjective elements in at least one of the stages. Thus, weighting method analysis shows that significant contributions can be made, since a large portion of these strategies still employ elements susceptible to error.
Only Shahraki and Noorossana [34] proposed a way of evaluating any variability parameter when selecting the best Pareto optimal solution. We therefore see the theoretical gap, that we intend to address in this study: the behavior of the prediction variance in relation to the weighting strategies. Now, consider the following problem [55]: where f i (x) are the objective functions to be optimized and w i are the weights assigned to each objective function.
In order to calculate the variance for the function described in Equation (20), the following process is considered: where ρ f i f j is the correlation between the functions f i and f j . Considering that it is possible to calculate the variance of f i (x) at a certain point (21) can be modified to: x 0 as constant for each function at a given point. Analyzing Equation (22), we see that the variance of the estimated responses is minimized by diversification, i.e., by the uniform distribution of weights among the functions involved in the MOP. Furthermore, negative correlations between responses tend to decrease variance. We propose the following: using entropic metrics to choose the optimal weights in MOP can reduce prediction variance. Now, consider the term as constant for each function at a given point. Analyzing Equation (22), we see that the variance of the estimated responses is minimized by diversification, i.e., by the uniform distribution of weights among the functions involved in the MOP. Furthermore, negative correlations between responses tend to decrease variance. We propose the following: using entropic metrics to choose the optimal weights in MOP can reduce prediction variance. Figure 1 shows the step-by-step proposal.
Step-by-step of the adopted methodology. Source: Adapted of Rocha et al. [55].

Experimental Design
The experimental data presented by Istadi and Amin [28] and Istadi and Amin [16] were used in this study. The authors sought to optimize a CO 2 -OCM process, by determining the condition that led to maximum CH 4 conversion, C 2 selectivity, and C 2 yield.
The CaO/CeO 2 and MnO/CeO 2 catalysts were prepared by impregnating ceria (CeO 2 ) with aqueous solutions of Ca(NO 3 ) 2 and Mn(NO 3 ) 2 , respectively, as described in Istadi and Amin [56] and Istadi and Amin [16]. The catalytic performances were tested in an experimental set-up as described in Istadi and Amin [56]. The CH 4 conversion, C 2 selectivity, and C 2 yield are calculated as defined by Wang et al. [26], considering that the carbon in methane is converted to C 2 H 6 , C 2 H 4 and CO.
The most frequently used experimental design for data collection, for modeling the response surface functions, is a Central Composite Design (CCD). According to Myers et al. [36], a CCD is chosen because it is an efficient design for sequential experimentation, allowing a reasonable amount of information to test the error without requiring a large number of experiments. Additionally, the design accommodates a spherical region with five levels for each factor, which is advantageous from an experimental region point of view. The CCD was then employed in the experimental design. Using CO 2 /CH 4 ratio (dimensionless), reactor temperature (K), wt.% CaO in ceria catalyst and wt.% MnO in ceria catalyst as the decision variables, a full factorial design 2 4 was performed with eight axial points and two center points, generating 26 experiments. The experimental matrix is shown in Table 1 in which x 1 is the CO 2 /CH 4 ratio, x 2 is the reactor temperature, x 3 is the wt.% CaO in ceria catalyst, x 4 is the wt.% MnO in ceria catalyst, y 1 is the CH 4 conversion (%), y 2 is the C 2 selectivity (%), and y 3 is the C 2 yield (%). The fixed experimental conditions are [16,28]: catalyst weight = 2 g; total feed flow rate = 100 mL/min; total pressure = 1 atm; gas hourly space velocity = 3000 mL/(g h). It is important to note that the experimental matrix presented in Table 1 was planned using only 2 CP. This can be very detrimental to the stability of the prediction variance. Myers et al. [36] do not recommend using a CCD with only two CP, because this practice does not guarantee good dispersion of the prediction variance throughout the experimental region, and analyzing how the proposed method behaves in such designs is important. This was why we chose this experimental matrix.
The decision variables were analyzed in coded form. They were decoded only at the end of the analyses. This was done using the following equation [51]: where: Hi and Lo are related to the values of level +1 and −1, respectively. The parameters used in the experiments and their levels are shown in Table 2. In the search for an optimization model that allows for the simultaneous conversion of CH 4 , C 2 selectivity , and C 2 yield, the decision-making process, based on multiple criteria, proposed by Rocha et al. [52] can be used. The authors sought to build a uniformly distributed Pareto frontier and a way to select the preferred Pareto optimal points as the ideal solution to a given problem. The authors used the NBI method, overcoming the disadvantages of the WSM method. To verify the robustness of the final result obtained, a variance metric was adopted. Although there are different measures of prediction performance to compare experimental designs, SPV is commonly adopted [57]. However, if direct comparisons between the expected variation of estimation are desired, UPV may be more appropriate.

Results and Discussion
The analysis of experimental data shown in Table 1 generated the mathematical modeling presented in Table 3. The values presented in bold represent the significant terms of the model (thus, p-value < 5%). Table 3 shows that the R 2 values indicate a good fit for the model. The lack-of-fit test is shown in Table 3. In the lack-of-fit test, small p-values are undesirable. At a 5% significance level, all presented models were adequate.
In order to establish a comparison of how each decision variable affects each response, the main effect plots for CH 4 conversion (y 1 ), C 2 selectivity (y 2 ) and C 2 yield (y 3 ) are shown in Figure 2.
According to this analysis, the reactor temperature (x 2 ) is the most significant factor in increasing CH 4 conversion (y 1 ). The reactor temperature (x 2 ) is also an important factor when analyzing C 2 selectivity (y 2 ) and C 2 yield (y 3 ), but in a different way. The reactor temperature (x 2 ) increases the C 2 selectivity (y 2 ) and C 2 yield (y 3 ) until reaching 1123 K. When increasing the temperature above 1123 K, C 2 selectivity (y 2 ) and C 2 yield (y 3 ) decrease. The drop is greater when analyzing the values of C 2 selectivity (y 2 ), and this could be attributed to more CH 4 being converted into CO rather than C 2 H 4 and/or C 2 H 6 [16]. For C 2 selectivity (y 2 ) and C 2 yield (y 3 ) the factors wt.% CaO in ceria catalyst (x 3 ) and wt.% MnO in ceria catalyst (x 4 ) showed very similar behavior, and the highest values for the responses obtained were extreme factor values, i.e., 5 and 25% CaO, and 1 and 9% MnO in ceria catalyst. The considerable impact that these factors have on C 2 selectivity (y 2 ) and C 2 yield (y 3 ) highlights the importance of the catalysts in promoting the product selectivity to C 2 H 4 and/or C 2 H 6, and in inhibiting the reaction to form CO and H 2 O [16]. The CO 2 /CH 4 ratio (x 1 ) had different behavior for C 2 selectivity (y 2 ) and C 2 yield (y 3 ). While the highest value for C 2 yield (y 3 ) was observed for the CO 2 /CH 4 ratio (x 1 ) 2, the highest values for the C 2 selectivity (y 2 ) obtained were in the extreme CO 2 /CH 4 ratio (x 1 ) range. The 2.5 CO 2 /CH 4 ratio (x 1 ) led to the highest CH 4 (y 1 ) conversion value with a drop in the value of this response at CO 2 /CH 4 ratio (x 1 ) 3. The abundance of CO 2 in the high CO 2 /CH 4 ratio most likely decreased the catalyst activity by covering the catalyst active sites [16]. When analyzing only reactor temperature (x 2 ), we saw that increased CH 4 conversion (y 1 ), C 2 selectivity (y 2 ) and C 2 yield (y 3 ) would be negatively affected, which corroborates the results presented by other authors [25][26][27]. We can see how these objectives are, thus, conflicting. The values presented in bold represent the significant terms of the model (thus, p-value < 5%). Table 3 shows that the R 2 values indicate a good fit for the model. The lack-of-fit test is shown in Table 3. In the lack-of-fit test, small p-values are undesirable. At a 5% significance level, all presented models were adequate.
In order to establish a comparison of how each decision variable affects each response, the main effect plots for CH4 conversion (y1), C2 selectivity (y2) and C2 yield (y3) are shown in Figure 2. According to this analysis, the reactor temperature (x2) is the most significant factor in increasing CH4 conversion (y1). The reactor temperature (x2) is also an important factor when analyzing C2 selectivity (y2) and C2 yield (y3), but in a different way. The reactor temperature (x2) increases the C2 selectivity (y2) and C2 yield (y3) until reaching 1123 K. When increasing the temperature above 1123 K, C2 selectivity (y2) and C2 yield (y3) decrease. The drop is greater when analyzing the values of C2 selectivity (y2), and this could be attributed to more CH4 being converted into CO rather than C2H4 and/or C2H6 [16]. For C2 selectivity (y2) and C2 yield (y3) the factors wt.% CaO in ceria catalyst (x3) and wt.% MnO in ceria catalyst (x4) showed very similar behavior, and the highest values for the responses obtained were extreme factor values, i.e., 5 and 25% CaO, and 1 and 9% MnO in ceria catalyst. The considerable impact that these factors have on C2 selectivity (y2) and C2 yield (y3) highlights the importance of the catalysts in promoting the product selectivity to C2H4 and/or C2H6, and in inhibiting the reaction to form CO and H2O [16]. The CO2/CH4  Figure 3 shows the response surfaces for CH 4 conversion (y 1 ), C 2 selectivity (y 2 ) and C 2 yield (y 3 ).
In order to check the convexity of the functions, before performing the multi-objective optimization, the nature of the stationary point was analyzed using Equation (14). For The eigenvalue signs indicate that the function is concave and the stationary point is a maximum point. Through an analysis of the nature of the stationary point, we see that the functions have different convexities. Thus, the weighted sum method for multi-objective optimization is not adequate [58]. We therefore adopted the NBI method. To implement the NBI optimization routine, initially, the payoff matrix was estimated, obtaining the results shown in Table 4. sponse at CO2/CH4 ratio (x1) 3. The abundance of CO2 in the high CO2/CH4 ratio most likely decreased the catalyst activity by covering the catalyst active sites [16]. When analyzing only reactor temperature (x2), we saw that increased CH4 conversion (y1), C2 selectivity (y2) and C2 yield (y3) would be negatively affected, which corroborates the results presented by other authors [25][26][27]. We can see how these objectives are, thus, conflicting. Figure 3 shows the response surfaces for CH4 conversion (y1), C2 selectivity (y2) and C2 yield (y3).    The payoff matrix is an underlying part of the NBI implementation. After this step, a mixture design for the weights of each objective function was defined, as presented in Table 5.   Table 5 presents the Pareto optimal set for the MOP under analysis. As shown in Rocha et al. [59], the prediction variance was affected by the weights assigned to the objectives. Table 6 shows the Pearson Correlation analysis between the values that were previously presented in Table 5. Using the Pearson Correlation analysis we see that the ξ weighting metric has a statistically negative and significant correlation with the UPV. Therefore, when maximizing this metric, the variance values tend to be lower. This information indicates that the search for the most preferred Pareto optimal point in multi-objective optimization using this metric leads to a robust response from a variability point of view. Using data presented in Table 5, the entropy modeling, GPE, weighting metric (ξ), and UPV were performed. Their canonical mixing polynomials (Equations (24)-(27)), response surfaces and contour graphics were obtained (Figures 4-7), which are presented below: UPV −0.228 0.642 −0.378 0.057 0.000 0.001 Using data presented in Table 5, the entropy modeling, GPE, weighting metric (ξ), and UPV were performed. Their canonical mixing polynomials (Equations (24)- (27)), response surfaces and contour graphics were obtained (Figures 4-7), which are presented below :   1  2  3  12  13   2 3  1 1 2 3  1 2 2 3  1 2 3 3   2  2  1 2  1  2  1 3  1  3  2    It is interesting to note that all canonical polynomials of mixtures showed a good fit, with R 2 close to 100%. The adopted variance metric, UPV, had the worst adjustments, with an R 2 value equal to 82.83%. However, this value may be acceptable [36]. This analysis highlights the fact that it is possible to model the variance metric in terms of weights. This is because the weights interfere in the solution space. However, since the optimization of distinct functions is performed simultaneously, the solution space is not the same as the initial area of the DOE, and therefore, its shape is distinguished from the Hat Matrix shape when modeling the variance.
Lastly, when maximizing ξ, as described in Equation (26), weights w 1 , w 2, and w 3 , related to the final solution, were identified as being w 1 = 0.2603, w 2 = 0.3202, and w 3 = 0.4195. These optimal weights were used in a multi-objective optimization of CH 4 conversion (y 1 ), C 2 selectivity (y 2 ), and C 2 yield (y 3 ), reaching 8.806%, 51.468%, and 3.275%, respectively. This result is different from that obtained by Istadi and Amin [28] and Istadi and Amin [16]. When optimizing the responses using a bio-objective optimization process, the authors neglected the real trade-off behavior between the responses. The lower diversification among the responses almost led to a single-response optimization. The greatest challenge of this process is to achieve both high CH 4 conversion (y 1 ) and high C 2 selectivity (y 2 ), since it has been proven that these responses are in conflict one with the other [25][26][27]. Despite this trade-off, the results are acceptable [25] especially with respect to using low-valued natural gas. Entropy 2021, 23, x FOR PEER REVIEW 16 of 22      It is interesting to note that all canonical polynomials of mixtures showed a good fit, with R 2 close to 100%. The adopted variance metric, UPV, had the worst adjustments, with an R 2 value equal to 82.83%. However, this value may be acceptable [36]. This analysis highlights the fact that it is possible to model the variance metric in terms of weights. This is because the weights interfere in the solution space. However, since the optimization of   It is interesting to note that all canonical polynomials of mixtures showed a good fit, with R 2 close to 100%. The adopted variance metric, UPV, had the worst adjustments, with an R 2 value equal to 82.83%. However, this value may be acceptable [36]. This analysis highlights the fact that it is possible to model the variance metric in terms of weights. This is because the weights interfere in the solution space. However, since the optimization of The optimal coded values of the decision variables are CO 2 /CH 4 ratio (x 1 ) = 1.0060, reactor temperature (x 2 ) = 0.7536, wt.% CaO in ceria catalyst (x 3 ) = 0.4358, and wt.% MnO in ceria catalyst (x 4 ) = 0.4822. Using Equation (23), the encoded values were transformed into uncoded values. Therefore, the optimal values of the decision variables are CO 2 /CH 4 ratio (x 1 ) = 2.50, reactor temperature (x 2 ) = 1179.5 K, wt.% CaO in ceria catalyst (x 3 ) = 17.2%, and wt.% MnO in ceria catalyst (x 4 ) = 6.0%. The high CO 2 /CH 4 ratio (x 1 ) favors CH 4 conversion (y 1 ) [26]. A high reactor temperature (x 2 ) leads to high CH 4 conversion (y 1 ), but the C 2 yield (y 3 ) decreases at reactor temperatures higher than 1173 K, due to low C 2 selectivity (y 2 ) [26,27]. The high wt.% CaO (x 3 ) and wt.% MnO (x 4 ) in ceria catalyst increases CH 4 conversion (y 1 ) and C 2 selectivity (y 2 ). According to Wang and Ohtsuka [27], these catalysts can activate CO 2 to produce active oxygen for CH 4 conversion (y 1 ) and their basicity leads to improvements in C 2 selectivity (y 2 ).
The data presented in Table 5 were used to construct Figure 8. Here, the Pareto frontier constructed using the NBI method is presented, with the ideal point highlighted.
The data presented in Table 5 were used to construct Figure 8. Here, the Pareto frontier constructed using the NBI method is presented, with the ideal point highlighted.  Figure 8 shows the uniform distribution of the optimal Pareto points on the frontier. This also presents the most preferrable Pareto point, which is the final solution for the MOP. The confidence interval values are shown in Table 7. For Optimal point Figure 8. Pareto frontier for CH 4 conversion (y 1 ), C 2 selectivity (y 2 ) and C 2 yield (y 3 ) optimization. Figure 8 shows the uniform distribution of the optimal Pareto points on the frontier. This also presents the most preferrable Pareto point, which is the final solution for the MOP. The confidence interval values are shown in Table 7. For α = 5%, the confidence interval of 100(1 − α)%, for the mean response, at point x 0 T = 1 x 01 x 02 . . . x 0k is: When analyzing the question involving variability, the final solution obtained that maximizes the ξ metric is robust, since this metric leads the solution in a region of minimum variance, with less variability, and greater reliability. However, the confidence interval of the C 2 selectivity (y 2 ) is higher than the other responses. This demonstrates that C 2 selectivity (y 2 ) is the most difficult parameter to control in the process. Figure 9 shows the overlap of the different objective functions defining the feasible region for the analyzed problem. Figure 9 shows the conflicting nature between CH 4 conversion (y 1 ) and C 2 selectivity (y 2 ). An increase in CH 4 conversion (y 1 ) leads to a decrease in C 2 selectivity (y 2 ). In this study, the optimum point was chosen based on the maximization of metric ξ, and still, it was shown to be the most robust point. However, one can see that the region of maximum C 2 yield (y 3 ) lies within the feasible region for the problem, and is obtainable, even though this is not the most reliable, from a statistical process control point of view.
When analyzing the question involving variability, the final solution obtained that maximizes the ξ metric is robust, since this metric leads the solution in a region of minimum variance, with less variability, and greater reliability. However, the confidence interval of the C2 selectivity (y2) is higher than the other responses. This demonstrates that C2 selectivity (y2) is the most difficult parameter to control in the process. Figure 9 shows the overlap of the different objective functions defining the feasible region for the analyzed problem. Figure 9. Point of optimization for CH4 conversion (y1), C2 selectivity (y2), and C2 yield (y3) optimization. Figure 9 shows the conflicting nature between CH4 conversion (y1) and C2 selectivity (y2). An increase in CH4 conversion (y1) leads to a decrease in C2 selectivity (y2). In this study, the optimum point was chosen based on the maximization of metric ξ, and still, it was shown to be the most robust point. However, one can see that the region of maximum C2 yield (y3) lies within the feasible region for the problem, and is obtainable, even though this is not the most reliable, from a statistical process control point of view.

Conclusions
This paper studied direct low-value natural gas conversion using carbon dioxide oxidative coupling of methane over CaO/CeO2 and MnO/CeO2 catalysts to produce C2 hydrocarbons. The NBI method was adopted to simultaneously optimize CH4 conversion, C2 selectivity, and C2 yield. This method helped build an evenly distributed Pareto frontier for the three responses, regardless of the convexity of the functions. Optimal point Figure 9. Point of optimization for CH 4 conversion (y 1 ), C 2 selectivity (y 2 ), and C 2 yield (y 3 ) optimization.

Conclusions
This paper studied direct low-value natural gas conversion using carbon dioxide oxidative coupling of methane over CaO/CeO 2 and MnO/CeO 2 catalysts to produce C 2 hydrocarbons. The NBI method was adopted to simultaneously optimize CH 4 conversion, C 2 selectivity, and C 2 yield. This method helped build an evenly distributed Pareto frontier for the three responses, regardless of the convexity of the functions. The mathematical response model showed acceptable fitting, and the functions were proven adequate. Some results described in literature have been corroborated, e.g., the positive influence of reactor temperature in CH 4 conversion, the negative influence of reactor temperature in C 2 selectivity, the synergistic effect between reactor temperature, the wt.% CaO in ceria catalyst in increasing CH 4 conversion and C 2 yield, and the synergistic effect between reactor temperature and wt.% MnO in ceria catalyst in increasing C 2 yield.
An entropic measure was used to select the most preferred ideal Pareto point for the final solution. The decision-making criteria was useful and fundamental in identifying and mapping regions with minimum variation within the optimal Pareto responses obtained in the optimization processes. Furthermore, this study demonstrates that the weights used in the multi-objective optimization process have an influence on the variation in the forecast of the responses obtained. This study also proves the robustness of the weighting process used in choosing the final solution.
The simultaneous optimal values for the objective functions were CH 4 conversion = 8.806%, C 2 selectivity = 51.468%, and C 2 yield = 3.275%. These results were obtained by using the following process parameter combinations: CO 2 /CH 4 ratio = 2.50, reactor temperature = 1179.5 K, wt.% CaO in ceria catalyst = 17.2%, and wt.% MnO in ceria catalyst = 6.0%. The analyses of the confidence interval for the responses showed that C 2 selectivity has greater variability, and was the most difficult parameter to control in the process.
From an environmental point of view, this is an efficient process since it helps reduce CO 2 and CH 4 emissions and reduces waste, as hydrocarbon resources, specifically from low-value natural gas, can be used.