Optimal Control of Hydrocarbons’ Hydrogenation with Catalysts

Abstract

In this paper, the optimal control problem of hydrocarbons’ hydrogenation was investigated in the presence of two catalysts—Nickel–Kieselguhr and Raney Nickel. This multistage chemical reaction holds significant practical importance, particularly in the production of high-density fuels. The optimal control problem was reformulated as a nonlinear global optimization problem and addressed using a modified Mind Evolutionary Computation algorithm. The proposed modifications include methods designed to ensure solution feasibility and ease of practical implementation. Using the proposed method, the performance of the two catalysts was compared under constant temperature conditions and with optimal control strategies. The results demonstrate that selecting an effective catalyst has a greater impact on the reaction’s efficiency than temperature control alone, with the Raney Nickel catalyst consistently outperforming the Nickel–Kieselguhr catalyst by at least 17%. Additionally, the optimization approach was applied to identify a new set of catalyst parameters. The newly obtained catalyst parameters allowed for the improvement of the results of the Raney Nickel catalyst by 18%. The results of all numerical experiments and implementation details are described in the paper.

1. Introduction

Many industrial processes, particularly in chemical production, require optimization to enhance efficiency, increase yields, and reduce the generation of undesirable by-products that may fail to meet environmental regulations. These complex, multistage processes are often modeled and simulated using systems of ordinary differential equations (ODEs), typically comprising tens of equations that can only be solved numerically [1,2,3].

The basis of the optimal control problem is an adequate mathematical model of the process kinetics. To solve the problem, it is important to develop effective computational algorithms and software packages for modeling using high-performance computing systems. Such systems are used in various fields, including modeling of two-phase mass transfer processes and other dissipative processes in fractured porous media [4].

One effective approach to addressing such challenges is through formulating an optimal control problem in which the desired behavior of a reaction is achieved by controlling one or more input parameters, such as reaction temperature or duration. These optimal control problems can be transformed into nonlinear global optimization tasks, which can then be solved using various optimization techniques, including evolutionary algorithms [5,6,7,8].

The application of evolutionary algorithms (EAs) to optimal control problems has gained significant traction in recent years due to their flexibility as well as their ability to quickly converge to suboptimal solutions, which are often sufficient in practice [9]. Many control problems involve nonlinear dynamics, time-varying constraints, and discontinuous behaviors. Evolutionary algorithms excel in these situations because they do not require gradient information and can handle complex, non-convex search spaces [10]. This is especially relevant for problems like spacecraft trajectory optimization, robotic control, or energy management, where system models are often too complex for analytical solutions [11]. Furthermore, EAs are well suited for optimizing multiple conflicting objectives, such as minimizing energy consumption while maximizing system performance—crucial for real-world applications like autonomous systems, where trade-offs between competing metrics are inevitable [12].

A key advantage of evolutionary methods is their adaptability, as they can incorporate problem-specific constraints, including time and resource limitations. Additionally, their ability to adapt dynamically to changes in the problem environment makes them robust for real-time or online optimization scenarios [10,11]. The parallelizable nature of EAs is another advantage, as it allows for simultaneous evaluation of multiple candidate solutions, making them highly efficient in modern computing environments [13,14].

Despite these strengths, evolutionary algorithms do have limitations, particularly when it comes to implementing solutions in practice. There is no guarantee that a control strategy derived from optimization will be feasible in a real-world setting. For instance, a solution might oscillate excessively or exhibit rapid changes, making it impractical for implementation [15].

In this paper, the problem of optimal control for a complex chemical reaction—the hydrogenation of hydrocarbons—is investigated in the presence of two catalysts: Nickel–Kieselguhr (NK) and Raney Nickel (RN). This process is of significant practical importance due to its role in jet fuel production and optimizing it could lead to considerable economic benefits [16]. The optimal control problem for this reaction was formulated and transformed into a global optimization task.

It should be noted that, due to the large number of available evolutionary algorithms, it can be challenging to make a reasonable selection of one algorithm over another. Even within a certain class of algorithms, the best algorithm is often problem-specific [15]. The Mind Evolutionary Computation (MEC) algorithm was selected as the primary optimization method due to its success in addressing similar problems in other studies [17]. MEC is particularly well suited for handling nonlinear optimization tasks, making it a fitting choice for this work [18]. In this study, a new modification of the traditional MEC algorithm is proposed to account for the practical aspects of control implementation discussed above. For instance, the obtained control function must align with real-world constraints, such as the discrete set of values a mixing device can handle or the feasible rate at which a reaction temperature can change [19]. The proposed approach ensures that such practical considerations are incorporated into the optimization process.

The main novelty of this study lies in the optimization algorithm proposed, which is based on the MEC algorithm for determining the optimal control of the catalytic hydrogenation process of polycyclic aromatic hydrocarbons, using a detailed kinetic model. The proposed method is applicable to various chemical processes when their kinetics are known, ensuring that the obtained solutions are feasible under given constraints.

The remainder of this paper is organized as follows: Section 2 presents the problem formulation, where the optimal control task is transformed into a nonlinear global optimization problem. Section 3 outlines the kinetic model of the hydrogenation process along with the corresponding optimization criteria. In Section 4, the modifications to the MEC algorithm are introduced, while Section 5 discusses the results of numerical experiments, analyzing them from computational and optimization perspectives. The conclusion summarizes the key findings and suggests directions for future research.

2. Problem Formulation

This paper deals with a global nonlinear constrained minimization problem.

$$\underset{X\in D\subset R^n}{\min }\Phi \left(X\right)=\Phi \left(X^{*}\right)={\Phi }^{*}$$

(1)

Here, $\Phi \left(X\right)$ is the scalar objective function, $X=\left(x_1,x_2,\dots ,x_n\right)$ is the n-dimensional vector of variables, $\Phi \left(X^{*}\right)={\Phi }^{*}$ is the required minimum value, $R^n$ is the n-dimensional arithmetical space, and D is the constrained search domain.

$$D=\left\{X|x_i^{min}\le x_i\le x_i^{max}, i\in \left[1:n\right]\right\}\subset R^n.$$

(2)

In our study, the control variable is the reaction temperature $T\left(t\right)$, where $t$ represents the reaction time. The goal is to find a control function $T\left(t\right)$ that optimizes the system’s behavior, such as maximizing the yield of target products. To convert this into a global optimization problem, the integration interval $\left[t_{start};t_{end}\right]$ was discretized into a set of discrete time points, with each interval $\left[t_i;t_{i+1}\right]$ corresponding to one second of reaction time. This time step is determined by the experimental unit, where temperature adjustments cannot occur more frequently.

The values of $T\left(t_i\right)$ represent the components of the vector of variables $X=\left(x_1,x_2,\dots ,x_n\right)$. In such a manner, the control function $T\left(t\right)$ is approximated with the use of a piecewise linear function that connects values $x_1,x_2,\dots ,x_n$. As a result, the objective function $\Phi \left(X\right)$ is transformed into the functional:

$$J\left(T\left(t\right)\right) \to \underset{T\left(t\right)}{\min }$$

(3)

which minimizes the difference between current and required behavior.

3. Hydrogenation of Hydrocarbons

Polycyclic aromatic hydrocarbons (PAHs) are a group of organic compounds characterized by molecules containing two or more benzene rings [20]. PAHs are found in heavy oil fractions [21,22] and can also be produced through laser irradiation of carbon materials. In practical applications, the presence of PAHs in hydrocarbon feedstocks negatively impacts performance characteristics. High-density fuels, for instance, require high density values (but not exceeding that of kerosene fractions) and a low concentration of aromatic hydrocarbons. These strict requirements contribute to the complexity and high cost of fuel production technologies [23,24,25].

The catalytic hydrogenation of PAHs, commonly referred to as dearomatization, was investigated in a laboratory setup using two nickel-based catalysts: Nickel–Kieselguhr (NK) and Raney Nickel (RN) [23,24]. This process achieved high yields, with up to $99\%$ of the product consisting of high-density components that meet jet fuel specifications. The resulting product exhibited key characteristics such as high density (up to $870 \mathrm{kg}/{\mathrm{m}}^3$), low aromatic hydrocarbon content, and the absence of sulfur. A detailed kinetic model for PAH hydrogenation in the presence of these catalysts was developed in [16]. Through kinetic experiments, the reaction model parameters for both catalysts were determined. The resulting kinetic model allows for conducting computational experiments to explore all possible reaction conditions, identifying the optimal ones. For each catalyst, a range of acceptable temperatures, reaction times, and temperature variations throughout the process can be evaluated numerically instead of through challenging and costly full-scale experiments.

The mathematical description of the change in concentrations of the process components was considered relative to the conditional time of contact of the reaction mixture with the catalyst. The equations of chemical kinetics, compiled to describe reactions, form a system of ordinary nonlinear differential equations in a closed system according to the law of mass action with the following initial conditions: $\tau =0, y_i\left(0\right)=y_i^0, \tau \in \left[0,{\tau }^{*}\right]$.

$${\displaystyle \frac{d{y}_i}{d\tau }}={\displaystyle \sum }_{j=1}^J{v}_{ij}{w}_j, i=1,\dots ,I,$$

(4)

$$w_j=k_j\left(A_j,E_j\right)·{\displaystyle \prod }_{i=1}^I{\left(y_i\right)}^{\left|{\alpha }_{ij}\right|}-k_{-j}\left(A_{-j},E_{-j}\right)·{\displaystyle \prod }_{i=1}^I{\left(y_i\right)}^{\left|{\beta }_{ij}\right|}$$

(5)

Here, $y_i$ is the concentration of the reaction reagent, $\mathrm{mole} \mathrm{fractions}; \tau$ is the conventional contact time of the reaction mixture with the catalyst, $\mathrm{kg}·\min /\mathrm{mol}; J$ is the number of stages; $I$ is the number of substances; $v_{ij}$ is the stoichiometric matrix; $w_j$ is the rate of the $j$-th stage, $1/\min$ or $\mathrm{mol}/\left(\mathrm{kg}·\min \right); k_j, k_{-j}$ are the rate constants of the stages, $1/\min$, determined by the Arrhenius equation; ${\alpha }_{ij}$ are the negative elements of the matrix; $v_{ij}, {\beta }_{ij}$ are the positive elements; $v_{ij}, A_j, A_{-j}$ are the pre-exponential factors, $1/\min$; $E_j, E_{-j}$ are the activation energies of the forward and reverse reactions, $\mathrm{kcal}/\mathrm{mol}$; and ${\tau }^{*}$ is the reaction duration, s.

Models (4) and (5) assume the Arrhenius equation for non-elementary stages. The main reactions of the process of the hydrogenation of monocyclic aromatic hydrocarbons and the hydrogenation of bicyclic aromatic hydrocarbons are detailed, including the identification of intermediate components. A detailed scheme of chemical transformations of the catalytic process of hydrogenation of polycyclic aromatic hydrocarbons is given in

Table 1. Scheme of chemical transformations and values of kinetic parameters of the catalytic process of hydrogenation of polycyclic aromatic hydrocarbons on a Nickel-on-Kieselguhr catalyst.

${\mathit{k}}_{\mathit{j}}$	Chemical Stages	${\mathit{E}}_{\mathit{j}}$	$\mathit{l}\mathit{n}\left({\mathit{A}}_{\mathit{j}}\right)$
$k_1$; $k_{21}$	$Y_1+3Y_2\leftrightharpoons Y_3$	15.42; 25.00	10.52; 6.00
$k_2$	$Y_1+Y_2\to Y_4+Y_5$	$6.75$	$3.68$
$k_3$; $k_{18}$	$Y_6+2Y_2\leftrightharpoons Y_7$	0.59; 16.30	0.45; 15.26
$k_4$	$Y_7+3Y_2\to Y_8$	$19.43$	$15.14$
$k_5$; $k_{19}$	$Y_9+3Y_2\leftrightharpoons Y_{10}$	10.18; 26.60	9.82; 24.23
$k_6$; $k_{20}$	$Y_{10}+3Y_2\leftrightharpoons Y_{11}$	31.60; 26.52	27.72; 24.33
$k_7$	$Y_{11}+Y_2\to Y_{12}+Y_{12}$	36.08	31.00
$k_8$	$Y_{10}+Y_2\to Y_4+Y_{12}$	25.99	23.04
$k_9$	$Y_4+3Y_2\to Y_{12}$	12.16	10.95
$k_{10}$	$Y_7+Y_2\to Y_{16}$	32.51	30.89
$k_{11}$	$Y_{16}+Y_2\to Y_1$	36.17	33.26
$k_{12}$	$Y_{16}+Y_2\to Y_{17}$	24.44	20.25
$k_{13}$	$Y_{12}\to Y_{18}$	28.54	25.59
$k_{14}$	$Y_{18}+Y_2\to Y_{19}$	34.42	32.16
$k_{15}$	$Y_{12}+Y_2\to Y_{20}$	37.18	32.42
$k_{16}$	$2Y_7\to Y_6+Y_8$	34.78	32.72
$k_{17}$	$2Y_9\to Y_9+Y_{11}$	32.48	28.53

.

In Table 1, $Y_1$—butylbenzene; $Y_2$—hydrogen; $Y_3$—butylcyclohexane; $Y_4$—benzene; $Y_5$—butane; $Y_6$—naphthalene; $Y_7$—tetralin; $Y_8$—decalin; $Y_9$—diphenyl; $Y_{10}$—cyclohexylbenzene; $Y_{11}$—dicyclohexyl; $Y_{12}$—cyclohexane; $Y_{16}$—methyllindane; $Y_{17}$—1-methyloctahydro-1-indene; $Y_{18}$—methylcyclopentane; $Y_{19}$—3-methylpentane; and $Y_{20}$—hexane. According to the scheme of chemical transformations (Table 1), the mathematical model of catalytic hydrogenation of polycyclic aromatic hydrocarbons is a system of 20 differential equations for changing the concentrations of substances according to Formulas (4) and (5) with 42 kinetic parameters.

To optimize the technological parameters for the hydrogenation of aromatic hydrocarbon concentrates, it is essential to understand how each parameter influences the conversion of raw materials. Kinetic factors, such as temperature and volumetric flow rate, form a distinct group of these parameters. The data gathered from these studies enabled the calculation of both the reaction kinetics and the required catalyst volume in the reactors. The impact of kinetic factors on the hydrogenation process was investigated for Nickel–Kieselguhr and Raney Nickel catalysts [16]. Using the developed kinetic model, it becomes necessary to formulate and solve the problem of optimal control for the catalytic PAH hydrogenation process.

In this reaction, there are two key target characteristics that need to be maximized, the conversion of the raw material $Z_1$ and the final concentration of naphthenes $Z_2$, at the end of the process. Since both criteria must be optimized simultaneously, the objective function takes the form of a linear convolution.

$$J\left(T\left(t\right)\right)=-\left({\mathrm{\lambda }}_1·Z_1+{\mathrm{\lambda }}_2·Z_2\right) \to \underset{T\left(t\right)}{\min }.$$

(6)

In this context, ${\mathrm{\lambda }}_1$ and ${\mathrm{\lambda }}_2$ are weighting coefficients used to balance the contributions of each criterion within the convolution. A negative sign is applied to convert the maximization problem into a minimization one. The system of differential equations consists of $20$ first-order equations, with the integration interval $t\in \left[0;120\right]$ in seconds. As previously mentioned, the time step for the entire reaction has been fixed at 1 s, yielding $121$ components in the variable vector $X$.

4. Modified Mind Evolutionary Computation Algorithm for an Optimal Control Problem

In this paper, the Mind Evolutionary Computation (MEC) algorithm [17] is examined. Since its introduction, the algorithm and its various adaptations [26,27] have been widely utilized in solving global optimization problems across multiple disciplines. MEC has proven effective in optimizing optical systems [18] and tackling optimal control issues in processes like the thermally stimulated luminescence of polyarylenphthalide [13] and the catalytic reforming of gasoline [28].

The original MEC algorithm models certain aspects of human behavior. Each individual $s$ is viewed as an intelligent agent operating within a group $S$ of similar individuals. Throughout the evolutionary process, every individual is influenced by others in the group. Achieving a prominent position within the group requires an individual to learn from the most successful members, while groups also follow a similar logic in intergroup competition.

The algorithm consists of three primary phases: group initialization, similar taxis, and dissimilation. The operations of similar taxis and dissimilation are repeated iteratively as the best value of the objective function $\Phi \left(X\right)$ improves. Once no further improvement is observed, the top-performing individual from the leading group is chosen as the solution to the optimization problem [29].

The population of the MEC algorithm consists of independent subpopulations with different instances of the SMEC algorithm. Each subpopulation is composed of leading groups ($S^b=\left(S_1^b,S_2^b,\dots ,S_{\left|S^b\right|}^b\right)$) and lagging groups ($S^w=\left(S_1^w,S_2^w,\dots ,S_{\left|S^w\right|}^w\right)$). The number of individuals within each group is set to be the same and equals $\left|S\right|$.

Every group has its own communication environment called a local blackboard, denoted as $C_k^l,$ $k\in \left[1:K\right]$, where $K$ is the number of groups. In addition, the whole multipopulation has a general global blackboard, i.e., $C^g.$ The global architecture of the MEC algorithm is presented in Figure 1. Here local and global blackboards contain vectors $X$, the corresponding values of the objective function $\mathrm{\Phi }\left(X\right)$ and their ranks $\#$ on the blackboard.

The canonical MEC algorithm has a few free parameters that can significantly affect its performance. In [30], the influence of the following free parameters was studied: the standard deviation ${\sigma }_{MEC}$ used for generating individuals; removal frequency of lagging groups ω, which is used to eliminate the most unsuccessful groups; and the ratio η between the number of lagging $S^w$ and leading groups $S^b$. In this study, the recommendations from [30] are used for the values of ω and η. For ${\sigma }_{MEC}$, a new adaptive strategy is proposed, which is discussed later.

In this work, several modifications are introduced to the canonical MEC algorithm to meet the practical demands of the reaction setup: a new step adaptation strategy, a constraint-handling mechanism utilizing the projection method [31], and a smoothing technique to generate feasible control strategies that are applicable in real-world scenarios. [28].

The description of the proposed modified algorithm is presented below.

• Initialization of groups within the search domain $D$ and problem decomposition.
  • Generate a specified number $\gamma$ of groups $S_i, i\in \left[1:\gamma \right]$; here, $\gamma$ is the free parameter of the algorithm.
  • Generate a random vector $X_{i,1}$. Identify this vector with the individual $s_{i,1}$ of the group $S_i$.
  • For every $X_i, i\in \left[1:\gamma \right]$, calculate the corresponding values of the objective function ${\Phi }_i$.
  • In each group, $S_i, i\in \left[1:\gamma \right]$, determine the initial coordinates of the rest of the individuals in group $s_{i,j}, j\in \left[2:\left|S\right|\right]$ according to the following formula:

    $$X_{k,i,j}=X_{k,i,1}+N_{\left|X\right|}\left(0,{\sigma }_{MEC}\right)$$

    (7)

Here, $\left|S\right|$ is the number of individuals in each group and ${\sigma }_{MEC}$ is the standard deviation used for initialization.

e.

Calculate the scores of all individuals in every group and put them on the corresponding local blackboards $C_i, i\in \left[1:\gamma \right]$.

f.

Create leading $S^b$ and lagging $S^w$ groups based on the obtained scores. The ratio between leading and lagging groups is determined by the free parameter $\eta .$ In this work, the number of lagging and leading groups is the same ($\eta =0.5$), corresponding to a situation when there is a balance between exploration and exploitation.

2.

The modified similar taxis stage.

• Take information on the current best individual $s_{i,j}, j\in \left[1:\left|S_i\right|\right]$ of the group $S_i$ from the blackboard $C_i$.
• Use the value of parameter ${\sigma }_{MEC}$ to generate new agent changes according to the following rule:

  $${\sigma }_{MEC}=\left\{\begin{array}{c}{\sigma }_0, if \lambda <\widehat{\lambda }, \\ {\displaystyle \frac{1}{{\left(\lambda -\widehat{\lambda }\right)}^{\theta }}}+\epsilon , if\lambda \ge \widehat{\lambda }\end{array}\right.$$

  (8)

Here, $\widehat{\lambda }$ is the threshold number of iterations, and when $\lambda \ge \widehat{\lambda }$, the standard deviation ${\sigma }_{MEC}$ begins to decrease; ${\sigma }_0$ is the initial value of the standard deviation; $\theta$ is the speed parameter (the recommended value $\theta =0.2$); and $\epsilon$ is the tolerance used to identify stagnation.

c.

Check for constraint violations:

• If any component of vector $X$ is outside of domain $D$, it is projected back onto the boundary of domain $D$;
• The difference between any two nearby components of vector $X$ is modified in order to not overcome the specified limit from physical experiments. This procedure is performed in random order either from $x_1$ to $x_n$ or vice versa. For this reaction, this limit is equal to $1 °\mathrm{C}$.

d.

Create new leading groups $S^b=\left(S_1^b,S_2^b,\dots ,S_{\left|S^b\right|}^b\right)$ and lagging groups $S^w=\left(S_1^w,S_2^w,\dots ,S_{\left|S^w\right|}^w\right)$ around current best individuals ${\tilde{s}}_{i,j}$ using Formula (7).

e.

Put information about the new winner on the corresponding local blackboard $C_k$ and the global blackboard $C_g$ available to the master.

3.

The dissimilation stage is performed after the similar taxis stage.

• Read the scores of all groups ${\Phi }_i^b, {\Phi }_j^w, i\in \left[1:\left|S^b\right|\right], j\in \left[1:\left|S^w\right|\right]$ from the global blackboard $C^g$;
• Compare all the scores. If the score of any leading group $S_i^b$ is less than the score of any lagging group $S_j^w$, then the lagging group becomes a leading group, and the first group becomes a lagging group. If the score of a lagging group $S_k^w$ is lower than the scores of all leading groups for $\omega$ consecutive iterations, then it is removed from the population.
• Replace each removed group with a new group via an initialization procedure.

For this study, the maximum allowed value of the objective function’s evaluation ${\lambda }_{max}$ was used as the termination criterion [30]. The block scheme of the proposed method is presented in Figure 2.

5. Computational Experiments and Analysis

In this study, we explored various aspects of the PAH hydrogenation process. First, we examined the influence of two different catalysts under constant temperature conditions without applying any control. Second, for both catalysts, we solved an optimal control problem for two criteria $Z_1$ and $Z_2$, both individually and as a linear convolution $J\left(t\right)$. The two catalysts exhibited significantly different behavior regarding the efficiency of the reaction. Finally, the proposed optimization approach was applied to the task of catalyst identification, aimed at discovering a new catalyst that could enhance the process’s output characteristics.

The modified MEC algorithm developed by the authors was implemented in Python. The following parameter values were used: the number of groups, $\gamma =20$; the number of individuals per group, $\left|S\right|=10$; the removing frequency, $\omega =30$ iterations; and the initial standard deviation, ${\sigma }_0=1$. For the experiments, a computational limit of $500$ iterations was set, resulting in $100,000$ objective function evaluations. All experiments were conducted on a virtual Windows server, utilizing a multi-start method with $20$ independent runs for each experiment.

As previously discussed, the MEC algorithm was selected for this study due to its proven effectiveness in related applications [14,17,18]. Benchmark studies [18,28] have compared the performance of the MEC algorithm, its various modifications, and other widely used optimization techniques, including the Particle Swarm Optimization (PSO) algorithm [31], the Gravity Search Algorithm (GSA) [32], and the Covariance Matrix Adaptation Evolution Strategy (CMAES) [33]. The primary innovation of this study lies in the proposed modifications that ensure the feasibility of the derived control solutions. A direct comparison of this approach with other methods that lack these modifications would be inherently biased, as it could result in solutions that, while potentially superior in terms of objective function values, are not practically implementable [13]. Conversely, incorporating these modifications into alternative algorithms would effectively result in a novel algorithmic framework, which, although potentially valuable for comparative analysis, extends beyond the scope of the current study.

5.1. Catalysts Comparison Under Constant Temperatures

The objective of this experiment is to assess the influence of the catalyst on the output characteristics, specifically the conversion of the raw material $Z_1$ and the final concentration of naphthenes $Z_2$, under constant reaction temperatures $T$. The feasible temperature range is $T\in \left[200;500\right] °\mathrm{C}$.

Figure 3 illustrates the change in raw material conversion $Z_1$ as a function of reaction temperature $T$ for two catalysts: Nickel–Kieselguhr (NK) and Raney Nickel (RN). Using the NK catalyst (blue curve), the maximum conversion is achieved at lower temperatures, between $230 °\mathrm{C}$ and $270 °\mathrm{C}$, with an exact value of $Z_1=0.709635$. At higher temperatures, conversion decreases slightly but remains above $0.65$. In contrast, the RN catalyst (red curve) reaches its peak performance at temperatures above $350 °\mathrm{C}$, with a maximum conversion of $Z_1=0.831832$, which is $17\%$ higher than the NK catalyst. Reaction temperature plays a more critical role with the RN catalyst, as conversion fluctuates between $0.5$ and $0.83$. In summary, the NK catalyst may be preferable when there are constraints on maximum reaction temperatures, particularly when $T<258 °\mathrm{C}$ (as shown by the intersection point in Figure 3).

Figure 4 presents similar results for the final concentration of naphthenes $Z_2$. The NK catalyst (blue curve) consistently yields lower concentrations of the target product across the entire range of reaction temperatures $T$. However, for both catalysts, the maximum concentration of naphthenes is achieved at the lowest temperature of $T=200 °\mathrm{C}$ ($Z_2^{NK}=0.023891$ and $Z_2^{RN}=0.066712$) Considering our goal is to maximize both $Z_1$ and $Z_2$, the Raney Nickel catalyst (red curve) appears more promising, as it produces higher concentration values at elevated temperatures. This suggests it can deliver a strong performance for both criteria.

5.2. Optimal Control with Different Catalysts

The second set of experiments was carried out to obtain an optimal control for the reaction with the use of the Nickel–Kieselguhr catalyst. First, we find an optimal control $T\left(t\right)$ for the conversion of the raw material $Z_1$ (Figure 5a). Utilizing such a control strategy instead of the constant temperature allows for increasing the value of $Z_1$ by roughly $1\%$ up to $Z_1=0.714042$. However, this control appears to be easy to implement in practice, as the temperature slowly decreases for the first 90 s of the reaction and then steadily increases up to $250 °\mathrm{C}$. This value is almost the same as the optimal constant temperature in the first set of experiments.

When dealing with the second criterion (Figure 5b), the final concentration of naphthenes $Z_2$, the obtained control represents an almost constant temperature except for the last 30 s of the reaction where the temperature $T\left(t\right)$ starts rising from $200 °\mathrm{C}$ to $205 °\mathrm{C}$. Such control allows for increasing the value of $Z_2$ by $0.5\%$ to $Z_2=0.02423$. Such control also appears to be practically feasible.

When optimizing the two criteria of $Z_1$ and $Z_2$ together as a linear convolution $J\left(t\right)$, the following weights were used: ${\mathrm{\lambda }}_1=1$ and ${\mathrm{\lambda }}_2=10$. Obtained controls are presented in Figure 6 in comparison with the individual controls. The corresponding values of the two criteria are $Z_1=0.686502$ and $Z_2=0.02423$. These results demonstrate that it is possible to simultaneously obtain high values for both criteria when using the Nickel–Kieselguhr catalyst. The value of $Z_1$ is slightly less than the best obtained value, while the value of $Z_2$ is the same.

The third set of experiments was conducted with the use of the Raney Nickel catalyst. In the same manner, we first study the conversion of the raw material $Z_1$ (Figure 7a). The obtained control $T\left(t\right)$ allowed for increasing the value of $Z_1$ also by $1\%$ up to $Z_1=0.837916$. The control strategy is also practically feasible as the temperature steadily increases up to $460 °\mathrm{C}$.

This temperature is significantly higher than the one required for the first catalyst. We decided to decrease the upper temperature limit for the RN catalyst to $300 °\mathrm{C}$. This would allow for the comparison of two catalysts when there are some energy restrictions. The yellow line in Figure 7a demonstrates that under specified conditions, it is necessary to maintain a constant temperature of $300 °\mathrm{C}$ in order to obtain the maximal value of $Z_1=0.782892$. This value is $10\%$ higher than the one obtained with the NK catalyst, thus proving the efficiency of the RN catalyst even under energy restrictions.

Regarding the second criteria $Z_2$ (Figure 7b), the obtained control is similar to the one obtained for the NK catalyst. The only difference is that a small temperature increase is required in the middle of the reaction. As a result, the output value of $Z_2$ was improved by $0.6\%$ up to $Z_2=0.067163$.

Obtained optimal control for two criteria as a convolution $J\left(t\right)$ also allows for providing high values for $Z_1$ and $Z_2$. For this experiment, we use the same weights of ${\mathrm{\lambda }}_1$ and ${\mathrm{\lambda }}_2$ as before. The results are presented in Figure 8 and demonstrate unexpected behavior. It is required to start the reaction at the highest possible temperature of $500 °\mathrm{C}$ and slowly decrease it after 40 s down to $450 °\mathrm{C}$. The resulting values of the two criteria are $Z_1=0.806898$ and $Z_2=0.039154$.

All obtained figures are summarized in

Table 2. Computational experiment results for investigated optimization criteria with a use of two catalysts.

Optimization Criterion	Nickel–Kieselguhr Catalyst	Raney Nickel Catalyst
Source Raw Material Conversion, $Z_1$	$Z_1=0.714042$; $Z_2=0.006386$;	$Z_1=0.837916$; $Z_2=0.019595$;
Target Concentration of Naphthenes, $Z_2$	$Z_1=0.68661$; $Z_2=0.02423$;	$Z_1=0.508969$; $Z_2=0.067163$;
Linear Convolution $J\left(t\right)$	$Z_1=0.686502$; $Z_2=0.02423$;	$Z_1=0.806898$; $Z_2=0.039154$;

. We can conclude that for the RN catalyst, two criteria are more sensitive than for the NK catalyst. The increase in $Z_2$ leads to a significant drop in $Z_1$ and vice versa. For the NK catalyst, the decrease in $Z_1$ is not that significant.

5.3. Catalyst Identification

The results from the previous sections clearly demonstrate that, for this reaction, selecting the right catalyst is more critical than controlling the reaction temperature $T$. Even with optimal control strategies, the NK catalyst cannot achieve the same output performance as the RN catalyst, emphasizing the practical importance of catalyst identification for process optimization.

In the model under investigation, a catalyst is characterized by 42 kinetic parameters. These parameters describe the chemical stages and influence the catalyst’s effect on the reaction. The activation energy $E_j$ defines the minimum energy required for the reaction, while the pre-exponential factor $A_j$ represents the total number of molecular collisions. Using these parameters (as shown in Table 1), we can calculate the rate constants for the stages based on the reaction temperature, following the Arrhenius equation.

We use these parameters as the input variable vector $X$ for our optimization approach, with the linear convolution $J\left(X\right)$ of $Z_1$ and $Z_2$ serving as the optimization criterion. The classical MEC algorithm for constrained global optimization problems is applied to address this problem, without incorporating the previously mentioned modifications for feasible control strategies.

The results are shown in Figure 9. In Figure 9a, we compare the RN catalyst’s parameters with the new optimized values obtained through the optimization process. Figure 9b highlights the differences between these parameters, indicating which should be increased and which should be decreased. The new set of parameters improved the value of $J\left(X\right)$, corresponding to the Raney Nickel catalyst, by 18%.

However, it is evident that the energy levels of the catalyst’s chemical components are not directly represented by the optimized values and require appropriate mapping to accurately identify the elements involved. This issue lies beyond the scope of the current study and will be addressed in future research.

6. Conclusions

In this study, the optimal control problem of hydrocarbons’ hydrogenation was addressed in the presence of two catalysts: Nickel–Kieselguhr and Raney Nickel. By reformulating the problem as a nonlinear global optimization task, the modified Mind Evolutionary Computation algorithm was proposed and utilized to optimize the reaction’s performance. The results clearly demonstrated that the choice of catalyst plays a more significant role in the efficiency of the reaction than controlling temperature alone. The RN catalyst consistently outperformed NK in both constant temperature conditions and optimal control scenarios, highlighting its superiority for practical applications in jet fuel production.

Moreover, by introducing modifications to the traditional MEC algorithm, it became possible to address practical constraints that often arise in real-world control problems, such as ensuring that the derived control functions are feasible for implementation. These modifications not only improved the robustness of the optimization but also provided a pathway for better integration of theoretical solutions into industrial setups.

Finally, the identification of new catalyst parameters using the proposed optimization approach resulted in significant performance improvements over the original RN catalyst. While the optimization outcomes were promising, additional research is needed to refine the mapping of these results to the actual energy levels of chemical substances in catalysts.

The insights gained from this study can be extended to other catalytic processes, offering a framework for optimizing complex, multistage chemical reactions in various industrial sectors. Further research will focus on refining the algorithm and exploring its application to a broader range of reactions and catalyst types.