Advancing Parameter Estimation in Differential Equations: A Hybrid Approach Integrating Levenberg–Marquardt and Luus–Jaakola Algorithms

Abstract

This study presents an integrated approach that combines the Levenberg–Marquardt (LM) and Luus–Jaakola (LJ) algorithms to enhance parameter estimation for various applications. The LM algorithm, known for its precision in solving non-linear least squares problems, is effectively paired with the LJ algorithm, a robust stochastic optimization method, to improve accuracy and computational efficiency. This hybrid LM-LJ methodology is demonstrated through several case studies, including the optimization of MESH equations in distillation processes, modeling controlled diffusion in biopolymer films, and analyzing heat and mass transfer during the drying of cylindrical quince slices. By overcoming the convergence issues typical of gradient-based methods and performing global searches without initial parameter bounds, this approach effectively handles complex models and closely aligns simulation results with experimental data. These capabilities highlight the versatility of this approach in engineering and environmental modeling, significantly enhancing parameter estimation in systems governed by differential equations.

1. Introduction

Optimization algorithms are essential tools in engineering and science for solving complex problems that require identifying the best set of parameters to minimize or maximize an objective function. Among the widely used algorithms in this context are the Luus–Jaakola (LJ) algorithm and the Levenberg–Marquardt (LM) algorithm.

The Luus–Jaakola algorithm is a direct search method known for its simplicity and effectiveness across a wide range of applications [1]. Its main advantage is that it does not require gradient information, which is particularly useful when computing derivatives is challenging or impractical [2]. The algorithm begins with a broad search to find an initial approximation of the global optimum, laying a solid foundation for refining the solution [3].

The LJ algorithm is versatile and has been applied in different fields, such as enhancing distributed generation in power systems [4], estimating aircraft parameters [5], scheduling irrigation systems [6], and optimizing control parameters for wing vibration systems [7]. It has also been enhanced by integrating elements from other algorithms to increase search efficiency and avoid premature convergence [8]. Moreover, it handles both continuous and discrete variables effectively, which is valuable for solving mixed-integer non-linear programming models that are common in real-world problems [9].

Despite its wide range of applications, the LJ algorithm has some limitations, especially in complex high-dimensional optimization problems, where exhaustive search can be slow and computationally intensive [3]. This challenge contrasts with hybrid methodologies, such as combining Luus–Jaakola with Hooke–Jeeves, which accelerate convergence by combining global search with local optimization strategies [10].

In contrast, the Levenberg–Marquardt (LM) algorithm is constructed for solving non-linear least squares problems by balancing between gradient descent and Gauss–Newton methods [11]. This balance is achieved through a damping parameter that adaptively selects the methodological approach based on the problem’s specific characteristics, thereby enhancing flexibility and convergence efficiency [12]. The algorithm effectiveness is greatly enhanced by using exact derivatives obtained through automatic differentiation. Combined with self-adapting regularization techniques and strategic limit tactics, these elements substantially accelerate the convergence of the LM algorithm [13].

This method is particularly advantageous in complex scenarios, such as optimizing backpropagation neural networks in non-linear systems, where it demonstrates superior processing accuracy and control effects compared to traditional methods [14]. Additionally, the LM algorithm has been adapted for a wide range of applications, from parameter identification in partial differential equations (PDEs) without needing the exact Fréchet derivative to tensor approximation problems [15]. It also effectively improves parameter calibration in inverter-based distributed energy resources, achieves high accuracy and rapid calibration [16], and enhances solution quality in heat conduction problems with singular scaling matrices [17,18].

Hybrid strategies combining the LJ algorithm have shown significant potential across various fields by enhancing optimization processes through the strengths of different methodologies. For example, in distributed generation systems, combining the Jaya algorithm with Luus–Jaakola (LJ) improves the convergence rate when selecting the location and size of distributed generations, considering demand-side management for system reliability and economic stability [2].

Similarly, hybrid strategies incorporating the LM algorithm for parameter estimation have demonstrated promise in various fields, showing improved performance, convergence speed, and accuracy. The LM algorithm, known for its robustness in solving non-linear least squares problems, has been effectively combined with other methods to address complex parameter estimation challenges. For instance, in Wiener–Volterra models, an optimization strategy using the LM algorithm is used to find kernel expansion parameters, achieving better results than previous works [19].

Both stochastic and deterministic methodologies have unique benefits that can be synergistically fused. One example is combining hybrid tactics in the LJ algorithm, which merges the broad exploration of stochastic techniques with the precise local exploitation provided by deterministic algorithms. This combination helps overcome the individual constraints of each method, such as the slow convergence of stochastic approaches and the risk of deterministic algorithms getting trapped in local optima.

This research aims to present a unified method that starts with the LJ algorithm for initial estimation, followed by refinement using the deterministic LM method. This approach is applied to chemical engineering problems governed by algebraic equations and systems of PDEs related to mass and heat transfer phenomena. This work demonstrates the integration of the LM and LJ algorithms, highlighting their effectiveness in parameter estimation across diverse applications. The efficiency of this combined approach is illustrated through various case studies, demonstrating significant improvements in parameter estimation accuracy and efficiency across different applications.

2. Methodology

2.1. Parameter Estimation

Least squares parameter estimation is employed to fit a model to data by minimizing the sum of squared differences between observed and predicted values. Various methods, including the forward–backward linear prediction total least squares (FB-TLS) method [20], the recursive least squares algorithm [21], and constrained minimization techniques [22], are applied for parameter estimation in diverse applications. These techniques are used to estimate parameters such as diffusion coefficients, system characteristics, and heat transfer coefficients, optimizing the model to best represent experimental data. The accuracy and effectiveness of these estimation techniques are validated through simulations and experimental studies, demonstrating improvements in spectral resolution, system identification, and model accuracy across different scenarios.

In this study, parameter estimation is performed when a mathematical model, containing unknown parameters, is supported by a dataset representing the system depicted by the model. Techniques for minimizing functions are employed to determine the values of these parameters, which often carry physical significance. The objective is to minimize the sum of squared errors, $S=\mathit{min}\left\{\sum _{i=1}^m\left(\sum _{j=1}^l{\left({\varphi }_{ij}^{experimental}-{\varphi }_{ij}^{model}\left(X{p}_n\right)\right)}^2\right)\right\}$, where $m$ is the number of equations; $l$ is the number of experimental data points; $n$ is the number of parameters to be estimated; $X{p}_n$ is the parameter to be estimated; ${\varphi }_{ij}^{model}$ is the data generated by the selected model; and ${\varphi }_{ij}^{experimental}$ is the experimental data.

2.2. Selection of Case Studies

This research explores five chemical and biochemical engineering case studies, including the hybrid strategy analysis in manufactured case studies. The investigation extends to an intricate analysis of MESH equations, which are crucial for designing and operating distillation systems [23]. Then, a pre-designed system of parabolic differential equations is analyzed. Furthermore, the study focuses on understanding diffusion phenomena within biopolymer films to improve material performance and quality [24]. Additionally, examining heat and mass transfer phenomena during the drying process of cylindrical quince slices is discussed, with an emphasis on identifying key variables that influence the efficiency and effectiveness of the drying operation [25].

2.3. Finite Difference Discretization of Non-Linear Partial Differential Equations

After identifying the model and the parameters to be estimated, the model is made dimensionless to facilitate spatial coordinate discretization using the finite difference method. This process transforms the partial differential equations (PDEs) into a system of ordinary differential equations (ODEs), where the number of ODEs corresponds to the number of nodes defined during discretization. Central finite differences are used for interior nodes, while forward and backward finite differences are applied to the boundary conditions representing the first and last nodes, respectively.

Figure 1 illustrates the implementation of finite difference discretization for a one-dimensional diffusion–reaction biological system. In this example, the function $\varphi$ depends on two independent variables and constants ($A_0, A_2,\dots , A_{10}$), which can serve as parameters ($X{p}_1,X{p}_2,\dots ,X{p}_{11}$). These parameters ($X{p}_n$), whether kinetic or transport-related, are defined according to the specific needs of the model representing the biological system.

2.4. Numerical Solution of Ordinary Differential Equation Systems

The numerical solution of ordinary differential equations in our model is performed using the fourth-order Runge–Kutta integration method. This technique applies numerical integration to the system equations ${\varphi }_{mod}\left(X_{pn}\right)$, with the step size $h$ adjusted according to the number of nodes used in the discretization. Specifically, when fewer than 15 nodes are involved, the recommended step size ranges from $1\times {10}^{-4}$ to $1\times {10}^{-3}$. For 15 or more nodes, a smaller step size, such as $1\times {10}^{-5}$ or smaller, is required to ensure accurate results.

2.5. Implementation of the Luus–Jaakola Algorithm

The Luus–Jaakola algorithm [1] is implemented to minimize the sum of squares represented by the function $S=min{\sum }_{i=1}^m\left({\sum }_{j=1}^l{\left({\varphi }_{ij}^{experimental}-{\varphi }_{ij}^{model}\left(X{p}_n\right)\right)}^2\right)$, where ${\varphi }_{ij}$ represents the experimental data, and ${\varphi }_{ij}^{model}\left(X{p}_n\right)$ is the numerical approximation of the model representing the system. Initial values for each parameter $(X{p}_1^0,X{p}_2^0,\dots ,X{p}_n^0)$ and an initial value for the search range $(r_1^0,r_2^0,\dots ,r_n^0)$ of each parameter are taken. A set of random values is generated for each parameter, denoted as $y_{ki}$. From this set, $p$ random values are selected, with $pn$ values for each of the $n$ parameters. A new value is assigned to each parameter according to the formula $X{p}_i^{j+1}=X{p}_i^j+y_{ki}{r}_i^j$, where $j=\mathrm{1,2},\dots ,j_{max}$, $i=\mathrm{1,2},\dots ,n$, and $l=\mathrm{1,2},\dots , p$.

Throughout each iteration, the algorithm generates new parameter values using a stochastic approach and adjusts the search range dynamically with a contraction factor $\left(ϵ\right)$. This iterative process continues until a pre-determined number of iterations is completed or the convergence criteria are met. The systematic reduction of the search range and continuous refinement of parameters underscore the rigor of the Luus–Jaakola algorithm, making it a robust approach for complex optimization problems. Once the algorithm completes the maximum number of iterations, the initial approximation of the parameter vector is obtained for use in the Levenberg–Marquardt algorithm. Algorithm 1 illustrates the pseudocode of the Luus–Jaakola algorithm.

Algorithm 1: Pseudocode for the implementation of the Luus–Jaakola algorithm

Start of the Luus–Jaakola Algorithm   Initialize parameters:     Define number of parameters/variables to estimate, $n$     Define objective function, $F\left(x\right)$           For solving system of $n$ equations, $F\left(x\right)=\mathit{min}\left\{\sum _{i=1}^n\left|f_i({Xp}_1,X{p}_2,\dots ,{Xp}_n)\right|\right\}$           For solving $n$ parameters, $F(x)=min{\sum }_{i=1}^m\left({\sum }_{j=1}^l{\left({\varphi }_{ji}^{experimental}-{\varphi }_{ij}^{model}\left(X{p}_n\right)\right)}^2\right)$     Define constrains, $g\left(x\right)\le 0$ and $h\left(x\right)\ge 0$     Define maximum number of iterations, $j_{max}$     Define maximum number of random values per iteration, $k_{max}$     Define contraction factor, $\epsilon$ (epsilon)     Define scaling factor, $f_{scaling}$     Initialize variable values $x_i^0=X{p}_i^0$ for each $i$ in $\{1, 2, \dots , n\}$     Initialize search range $r_i^0$ for each $i$ in $\{1, 2, \dots , n\}$   For $j=1$ to $j_{max}$ do:     For $i=1$ to $n$ do:       For l $=1$ to $p$ do:         Generate random value, $Y_{i,l}$, within the range $[x_i^j-r_i^j, x_i^j+r_i^j]$         Evaluate constraints $g_i\left(Y_{i,l}\right)\le 0$ and $h_j\left(Y_{i,l}\right)\ge 0$         If all constraints are satisfied:           Evaluate the objective function $F\left(x\right)$ at $Y_{i,l}$           If $F\left(Y_{i,l}\right)$ improves the current solution:             Update the best solution found with $Y_{i,l}$     Update the initial set for the next iteration:       For $i=1$ to $n$ do:         $x_i^{j+1}$ = Best value found for $x_i^j$ that minimized $F\left(x\right)$         Reduce search range: $r_i^{j+1}=(1-\epsilon )\mathrm{\hspace{1em} }{r}_i^j$         Save $x_i^{best}=x_i^{j+1}$   Display the best solution found, $x_i^{best}$ that minimized $F\left(x\right)$ End

2.6. Implementation of the Levenberg–Marquardt Algorithm

Using the initial vector approximation provided by the Luus–Jaakola algorithm, the numerical approximation of the parameter vector is determined using the deterministic Levenberg–Marquardt method. The Levenberg–Marquardt method is defined as follows: $[X_{pn}^{\left(i+1\right)}=X_{pn}^{\left(i\right)}-{\left[H_m+\lambda \hspace{0.17em}\mathit{diag}\left(H_m\right)\right]}^{-1}\nabla S]$, where $(\nabla S)$ is the gradient of the function $S=min {\sum }_{i=1}^m{\sum }_{j=1}^l{\left({\varphi }_{\mathrm{ij}}^{experimental}-{\varphi }_{\mathrm{ij}}^{model}\left(X_{pn}\right)\right)}^2$ with respect to each parameter to be estimated ($X_{pn}$). The gradient $\nabla S$ is represented as $\left\{{\displaystyle \frac{\partial S}{\partial X_{p1}}},{\displaystyle \frac{\partial S}{\partial X_{p2}}},\dots ,{\displaystyle \frac{\partial S}{\partial X_{pn}}}\right\}$. The Hessian matrix $\left(H_m\right)$ is an $( n\times n )$ matrix, where $( n )$ corresponds to the number of parameters to be estimated and is calculated as $H_m=\left[{\displaystyle \frac{{\partial }^{2}S}{\partial X_{pi}\partial X_{pj}}}\right]$. The term $\lambda \cdot \mathit{diag}\left(H_m\right)$ refers to the Levenberg–Marquardt parameter $\lambda$, which scales the diagonal elements of the Hessian matrix, defined as $\lambda \cdot \mathit{diag}\left(H_m\right)=\left[\lambda {\displaystyle \frac{{\partial }^{2}S}{\partial X_{p1}^2}},0,\dots ,0;0,\lambda {\displaystyle \frac{{\partial }^{2}S}{\partial X_{p2}^2}},\dots ,0\right]$. Algorithm 2 illustrates the pseudocode of the Levenberg–Marquardt algorithm [11].

Algorithm 2: Pseudocode for the implementation of the Levenberg–Marquardt algorithm

Start of the Levenberg–Marquardt Algorithm   Initialize parameters:     Set initial LAMBDA $\lambda =1$     Define convergence tolerance, $tolerance=1\times {10}^{-5}$     Set maximum number of iterations, $ma{x}_{iterations}$     Initialize the parameter vector, $X{p}_n^i=x_i^{best},$ with values from Luus–Jaakola that minimized $F\left(x\right)$             For solving system of $n$ equations, $F\left(x\right)=\mathit{min}\left\{\sum _{i=1}^n\left|f_i({Xp}_1,{Xp}_2,\dots ,{Xp}_n)\right|\right\}$             For solving parameter estimation, $F(x)=min{\sum }_{i=1}^m\left({\sum }_{j=1}^l{\left({\varphi }_{ij}^{experimental}-{\varphi }_{ij}^{model}\left(X{p}_n\right)\right)}^2\right)$   For $i=1$ to $ma{x}_{iterations}$ do:     Calculate the gradient of $F$, $\nabla F={\left[{\displaystyle \frac{\partial F}{\partial X{p}_1}}, {\displaystyle \frac{\partial F}{\partial X{p}_2}},\dots ,{\displaystyle \frac{\partial F}{\partial X{p}_n}}\right]}^T$     Calculate the Hessian matrix, $H_m=\left[\begin{array}{ccc}{\displaystyle \frac{{\partial }^{2}F}{\partial X{p}_1^2}}& \begin{array}{cc}{\displaystyle \frac{{\partial }^{2}F}{\partial X{p}_1\partial X{p}_2}}& \dots \end{array}& {\displaystyle \frac{{\partial }^{2}F}{\partial X{p}_1\partial X{p}_n}}\\ {\displaystyle \frac{{\partial }^{2}F}{\partial X{p}_2\partial X{p}_1}}& {\displaystyle \frac{{\partial }^{2}F}{\partial X{p}_2^2}}& {\displaystyle \frac{{\partial }^{2}F}{\partial X{p}_2\partial X{p}_n}}\\ \begin{array}{c}⋮\\ {\displaystyle \frac{{\partial }^{2}F}{\partial X{p}_n\partial X{p}_1}}\end{array}& \begin{array}{cc}\begin{array}{c}⋮\\ {\displaystyle \frac{{\partial }^{2}F}{\partial X{p}_n\partial X{p}_2}}\end{array}& \begin{array}{c}\ddots \\ \cdots \end{array}\end{array}& \begin{array}{c}⋮\\ {\displaystyle \frac{{\partial }^{2}F}{\partial X{p}_n^2}}\end{array}\end{array}\right]$     Adjust the $H_m$ with $\lambda$: $H_m=H_m+\lambda diag\left(H_m\right)$     Solve the equation system $H_m \ast b=\nabla S$ using LU factorization.     Update the parameter vector: $X{p}_n^{i+1}=X{p}_n^i-{\displaystyle \frac{\nabla S}{[H_m+\lambda diag(H_m\left)\right]}}$     Calculate $F$ using the new parameters, $X{p}_n^{i+1}$     If $F^{i+1}<F^i$ then:       Reduce λ by a factor, e.g., 0.9     Else:       Increase λ by a factor, e.g., 1.1     Check for convergence:       If $abs(\nabla F) \le tolerance$ for all components then:         Stop iterations.     Display the current parameter vector, $X{p}_n^{i+1}$ End

2.7. Computation Algorithm

The Luus–Jaakola and Levenberg–Marquardt computational algorithms were developed in FORTRAN^® using Intel 64 Visual Studio 2019 environment on a computer with an Intel (R) Core (TM) i9–7900X CPU @ 3.30 GHz, 32 GB RAM, and a 64-bit operating system and x64 processor (AZTX SISTEMAS, Celaya, Mexico). To generate visual representations of the operating conditions of the algorithms within the hybrid strategy, the hybrid algorithms were developed in MATLAB^® R2023b on a computer with an Intel (R) Core (TM) i7-9750H CPU @ 2.60GHz, 16.0 GB of RAM, and a 64-bit Windows 10. The structure of the algorithms developed in MATLAB^® and FORTRAN^® follows the flowchart shown in Figure 2.

3. Results and Discussion

The hybrid optimization strategy, integrating the Luus–Jaakola (LJ) and Levenberg–Marquardt (LM) algorithms, has proven effective across a range of chemical engineering applications, as demonstrated in the following case studies. The methodology effectively addresses challenges related to initial approximation and convergence, which are common in complex optimization problems, particularly those involving non-linear systems and differential equations.

3.1. Case Study 1: Solution of MESH Equations in Distillation Column

A simple distillation column with five stages is described for separating propane, n-butane, and iso-butane. The feed stream (stage 3) is ${\mathrm{z}}_{(\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{a}\mathrm{n}\mathrm{e},3)}=0.10$, ${\mathrm{z}}_{(\mathrm{n}-\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{e},3)}=0.70$, and ${\mathrm{z}}_{(\mathrm{i}\mathrm{s}\mathrm{o}-\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{e},3)}=0.20$. The boiler provides 2691.27 J/s of heat, with the distillate extracted at 0.0001633 kgmol/h and the bottom product at 0.0002903 kgmol/h. The vapor–liquid equilibrium ratios for component $i$ at stage $j$ are $K_{\left(i,j\right)}=T_i{\left(a_{1,i}+a_{2,i} T_j+a_{3,i}{T}_j^2+a_{4,i}{T}_j^3\right)}^3$, where the values of the constants are taken from [23].

Consider the liquid–vapor equilibrium equations (E-equations) $y_{i,j}=K_{i,j}{x}_{i,j}$; the material balance is given by $L_{j-1}{x}_{i,j-1}+F_j{z}_{i,j}+V_{j+1}{y}_{i,j+1}=\left(L_j+U_j\right)x_{i,j}$; the enthalpy balance is given by $Q+L_{j-1}{h}_{j-1}+F_j{h}_{F,j}+V_{j+1}{H}_{j+1}=\left(L_j+U_j\right)h_j$, with the sum of molar fractions (S-equations) $\sum _{i=1}^{n_c}{y}_{i,j}, y \sum _{i=1}^{n_c}{x}_{i,j}$, where one of the S-equations is taken from the total mole balance equation $L_{j-1}+F_j+V_{j+1}=L_j+U_j+V_j+W_j$.

With the MESH equations defined, an algorithm is developed to compute the 27 equations that describe the distillation problem. These equations were used for optimization with two different tools: in FORTRAN, optimization was performed using the Luus–Jaakola algorithm, and in MATLAB, the ‘fsolve’ function was used by applying the Levenberg–Marquardt algorithm. In this case, it was identified that there was difficulty in achieving convergence in the system of equations due to the initial approximation used.

In the case of the Levenberg–Marquardt algorithm, the initial approximation must necessarily be close to the solution of the system to find convergence; with initial approximations of 0, 1, 2, and 3, no solution to the system is found. However, when using the Luus–Jaakola algorithm starting from approximations of 0, 1, 2, and 3, an approximation is found that can be refined with the Levenberg–Marquardt method and converges to the solution of the system.

It is worth highlighting that the primary use of the Luus–Jaakola algorithm in complex engineering problems is to find an approximation that can be refined with a deterministic method to find the solution of the study system. It is crucial to note that the main challenge in these case studies is finding an optimal initial approximation when there is a wide range of independent variables used. In this instance, it is notable that the temperature values range from 230 to 650 K, while the compositions of liquid and vapor, $\left(x_{\left\{i,j\right\}}\right)$ and $\left(y_{\left\{i,j\right\}}\right)$, are fractions less than 1.

When the variability in magnitudes complicates the establishment of an initial approximation that illustrates smooth convergence, it is imperative to include an extra scaling parameter in the Luus–Jaakola algorithm, given a broad comprehension of the variable magnitudes to be approximated. This adjustment guarantees that the estimation of all variables is aligned with $X{p}_i^j=\left(X{p}_i^{j-1}+y_{ki}{r}_i^{j-1}\right)\ast \left(f_{scaling}\right)$, wherein the scaling factor is modified to correspond with the magnitudes of the variables present in the MESH equations.

Table 1. Results for the initial and final stages of the distillation column using Luus–Jaakola/Levenberg–Marquardt strategy.

Algorithm	$\mathbf{Solution} {\left[{\mathbf{x}}_1^{\mathbf{s}\mathbf{o}\mathbf{l}},\dots ,{\mathbf{x}}_{\mathbf{n}}^{\mathbf{s}\mathbf{o}\mathbf{l}}\right]}^{\mathbf{T}}$	${\displaystyle \sum _{\mathbf{i}=1}^{\mathbf{n}}}\left|\mathbf{f}\left({\mathbf{x}}_1^{\mathbf{s}\mathbf{o}\mathbf{l}},\dots ,{\mathbf{x}}_{\mathbf{n}}^{\mathbf{s}\mathbf{o}\mathbf{l}}\right)\right|$	CPU Time
Luus–Jaakola	Stage 1 [${\mathrm{x}}_{\left(\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{a}\mathrm{n}\mathrm{e},1\right)}=0.27060,{\mathrm{x}}_{\left(\mathrm{n}-\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{e},1\right)}=0.28049,{\mathrm{x}}_{\left(\mathrm{i}\mathrm{s}\mathrm{o}-\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{e},1\right)}=0.25871,{\mathrm{T}}_1=550.74751\left(\mathrm{K}\right)]$ … Stage 5 $[{\mathrm{x}}_{\left(\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{a}\mathrm{n}\mathrm{e},5\right)}=0.26826,{\mathrm{x}}_{\left(\mathrm{n}-\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{e},5\right)}=0.29115,$ ${\mathrm{x}}_{\left(\mathrm{i}\mathrm{s}\mathrm{o}-\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{e},5\right)}=0.27540, {\mathrm{T}}_5=551.06654\left(\mathrm{K}\right)]$	2906.47366	9 s
Levenberg–Marquardt	Without convergence	Without convergence	-----
Hybrid strategy	Stage 1 [${\mathrm{x}}_{\left(\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{a}\mathrm{n}\mathrm{e},1\right)}=0.1062,{\mathrm{x}}_{\left(\mathrm{n}-\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{e},1\right)}=0.6847,{\mathrm{x}}_{\left(\mathrm{i}\mathrm{s}\mathrm{o}-\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{e},1\right)}=0.2091,$ ${\mathrm{T}}_1=680.8999\left(\mathrm{K}\right)]$ … Stage 5 [${\mathrm{x}}_{\left(\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{a}\mathrm{n}\mathrm{e},5\right)}=0.0429,{\mathrm{x}}_{\left(\mathrm{n}-\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{e},5\right)}=0.7685, {\mathrm{x}}_{\left(\mathrm{i}\mathrm{s}\mathrm{o}-\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{e},5\right)}=0.1885,{\mathrm{T}}_5=690.1339\left(\mathrm{K}\right)]$	2.67 × 10−12	51.6 s

Luus–Jaakola algorithm: 20,000 maximum iterations, 2000 random values per iteration. Levenberg–Marquardt algorithm: 30,156 maximum iterations, 859,475 evaluations of $\mathrm{f}\left({\mathrm{x}}_1,{\mathrm{x}}_2,\dots ,{\mathrm{x}}_{\mathrm{n}}\right)$. Operating conditions of the hybrid strategy. Luus–Jaakola algorithm: 20,000 maximum iterations, 2000 random values per iteration. Levenberg–Marquardt algorithm: 10,319 maximum iterations, 294,115 evaluations of $\mathrm{f}\left({\mathrm{x}}_1,{\mathrm{x}}_2,\dots ,{\mathrm{x}}_{\mathrm{n}}\right)$.

presents the results for both the initial and final stages of the distillation column, illustrating the effectiveness of the Luus–Jaakola (LJ) algorithm, the Levenberg–Marquardt (LM) method, and their hybrid strategy in parameter estimation. The LJ algorithm provided initial estimates with a sum of absolute errors of 2906.47366 and a CPU time of 9 s. In contrast, the LM algorithm failed to converge when used independently. The hybrid strategy, however, significantly improved the parameter estimates, reducing the sum of absolute errors to $2.67\times {10}^{-12}$ and achieving a highly accurate solution. However, this enhanced precision came with an increased CPU time of 51.6 s. These results underscore the hybrid strategy’s ability to handle complex optimization problems, effectively manage convergence issues, and deliver highly accurate parameter estimations. Figure 3 illustrates the application of the LJ algorithm in resolving MESH equations, highlighting the results obtained at various iteration points.

To analyze the performance of the direct search executed with the Luus–Jaakola algorithm, various operating conditions were examined to evaluate the possibility of optimizing the algorithm. The operating conditions evaluated were the size of the contraction factor (ϵ), the search range, the number of iterations, and the number of random values generated.

Figure 3 illustrates the application of the LJ algorithm in resolving MESH equations, resulting in a thorough analysis of the outcomes achieved at various iteration points, such as the 1st, 15th, and 50th iterations. These results cover a range of crucial elements that contribute to the comprehensive analysis of the Luus–Jaakola algorithm. Initially, they incorporate the measurement of the molar ratio of liquid propane at the beginning of the process, providing insight into the essential initial conditions. Subsequently, they include the assessment of the exit temperature, a critical factor impacting the overall effectiveness and operation of the system. Finally, the results incorporate the accurate determination of the flow rate at the start of the system, which is essential for understanding the dynamics and behavior of the process.

As the iterations advance within the Luus–Jaakola algorithm, the search range gradually diminishes as it approaches the minimum of the function. Consequently, the values produced in each iteration refine the function of seeking a local minimum at the very least. As illustrated in Figure 3, increasing the number of iterations results in a lower minimum value of the objective function. This occurs as new random values are generated, and the search range decreases.

The scaling factors were generated to narrow down the system solution, and it is essential to modify them when there is a disparity in the order of magnitude of the system variables. This ensures better approximation to the system solution using the LJ algorithm. This case study highlights the potential of using the Luus–Jaakola algorithm in distillation problems, where temperatures $\left(T\right)$, liquid compositions $\left(x_{i,j}\right)$ and vapor compositions $\left(y_{i,j}\right)$, and flows on each tray are estimated based on the generation of MESH equations. This case study emphasizes the importance of using the LJ algorithm to generate an initial approximation that can be refined with the LM algorithm, especially in complex engineering problems like distillation, where the initial approximation must be close to the solution for successful convergence.

3.2. Case Study 2: Estimation of Four Parameters in a Manufactured PDE Set

The manufactured case is composed of three analytical solutions (Equations (1)–(3)).

$$T_1=X_{p1}\mathit{sin}\left(x\right)+{\displaystyle \frac{2x}{t+X_{p2}}}$$

(1)

$$T_2=X_{p3}x e^{-t}$$

(2)

$$T_3=x^3+{\displaystyle \frac{2x}{t+X_{p4}}}$$

(3)

From the analytical solution, the parameters $X{p}_1,X{p}_2,X{p}_3$, and $X{p}_4$ are obtained with a solution vector $X{p}_n^T={\left[\mathrm{3.0,1.0,2.0,0.5}\right]}^T$. A system of three partial differential parabolic equations is designed, each with its corresponding boundary conditions and initial condition with a different structure per differential equation. The first PDE is described in Equation (4).

$${\displaystyle \frac{\partial T_1}{\partial t}}=2\left({\displaystyle \frac{{\partial }^2{T}_1}{\partial x^2}}\right)+{\left({\displaystyle \frac{\partial T_2}{\partial x}}\right)}^2-{\displaystyle \frac{2x}{{\left(t+X_{p2}\right)}^2}}+2X_{p1}\mathit{sin}\left(x\right)-X_{p3}{e}^{-t}$$

(4)

where the initial condition is set at $t=0,T_1=X{p}_{1}sin\left(x\right)+2x/X{p}_2$, with boundary condition 1 at $x=0, T_1=0$ and boundary condition 2 at $x=1, T_1=3sin\left(1\right)+2x/\left(t+X{p}_2\right).$ The second PDE is described in Equation (5).

$${\displaystyle \frac{\partial T_2}{\partial t}}=3\left({\displaystyle \frac{{\partial }^2{T}_2}{\partial x^2}}\right)-X_{p3}x{e}^{-t}$$

(5)

where the initial condition is set at $t=0,T_2=X{p}_{3}x$, with boundary condition 1 at $x=0,T_2=0$ and boundary condition 2 at $x=1,T_2=2e^{\left(-t\right)}.$ Finally, Equation (6) describes the third PDE for this manufactured case.

$${\displaystyle \frac{\partial T_3}{\partial t}}={\displaystyle \frac{\partial T_1}{\partial x}}+{\displaystyle \frac{{\partial }^2{T}_3}{\partial x^2}}-{\displaystyle \frac{2x}{{\left(t+X_{p4}\right)}^2}}-X_{p1}\mathit{cos}\left(x\right)-{\displaystyle \frac{2}{t+X_{p2}}}-6x$$

(6)

where the initial condition is set at $t=0,T_3=x^3+2x/X{p}_2$, with boundary condition 1 at $x=0,T_3=0$ and boundary condition 2 at $x=1,T_3=1+2/\left(t+X{p}_4\right)$.

The experimental data generated for parameter estimation are obtained from the analytical solution at $x=0.5,X{p}_n^T={\left[3.0, 1.0, 2.0, 2.0, 0.5\right]}^T$, over a time range from 0 to 1. For handling the system of partial differential equations (PDEs), the spatial domain of each PDE was discretized along with their corresponding boundary conditions, resulting in three sets of ordinary differential equations (ODEs). These systems of equations were structured into mesh configurations with 7, 15, 21, and 35 nodes, respectively. The solutions to these systems were obtained using the fourth-order Runge–Kutta method. The process involves adjusting the shrinkage factor, particularly when the number of nodes varies, during parameter estimation using the Luus–Jaakola algorithm.

Table 2. Parameter estimation using the Luus–Jaakola algorithm for mesh analysis of a manufactured ODE set.

Nodes	Solution	Quadratic Sum of Residues (S)	CPU Time (s)
7	[2.25703, 0.61139, 1.16794, 0.51272] a1	3.08889	3.9
	[2.88156, 0.90639, 1.99978, 0.49958] a2	0.01843	3.9
	[2.99828, 0.99928, 1.99524, 0.49996] a3	0.00004	3.9
	[2.25702, 0.61139, 1.16792, 0.51273] b1	3.08903	2.0
	[2.87198, 0.89723, 2.00609, 0.49876] b2	0.02246	2.0
	[2.96524, 0.97154, 1.98258, 0.50055] b3	0.00218	2.0
15	[2.25796, 0.61267, 1.16962, 0.51642] a1	3.04299	5.5
	[2.85905, 0.89233, 1.98194, 0.50124] a2	0.02591	5.5
	[2.99448, 0.99971, 1.98668, 0.50146] a3	0.00027	5.5
	[2.25795, 0.61267, 1.16960, 0.51643] b1	3.04312	2.8
	[2.84861, 0.88345, 1.98602, 0.49980] b2	0.03068	2.8
	[2.96929, 0.98880, 1.97909, 0.50283] b3	0.00223	2.8
21	[2.25801, 0.61280, 1.16964, 0.51659] a1	2.98875	205.4
	[2.93414, 0.95307, 1.97235, 0.50324] a2	0.00555	238.7
	[3.00096, 1.00956, 1.97129, 0.50273] a3	0.00130	138.1
	[2.25800, 0.61281, 1.16962, 0.51660] b1	2.98888	67.3
	[2.92489, 0.94595, 1.97954, 0.50200] b2	0.00737	69.1
	[2.98567, 0.98985, 1.97134, 0.50100] b3	0.00200	69.1
35	[2.25247, 0.61618, 1.17278, 0.52120] a1	2.78752	193.7
	[2.96340, 0.99266, 1.91425, 0.50918] a2	0.01140	192.4
	[2.98273, 1.00899, 1.91535, 0.50900] a3	0.01134	196.5
	[2.25246, 0.61619, 1.17277, 0.52121] b1	2.78764	94.5
	[2.95631, 0.98824, 1.91714, 0.50943] b2	0.01164	122.6
	[2.93932, 0.97236, 1.92994, 0.50622] b3	0.01417	92.6

Operating conditions of the Luus–Jaakola algorithm: ^a1 100 maximum iterations, 20 random values per iteration, $h=1\times {10}^{-4}$ and є = 0.80. ^a2 100 maximum iterations, 20 random values per iteration, $h=1\times {10}^{-4}$ and є = 0.90. ^a3 100 maximum iterations, 20 random values per iteration, $h=1\times {10}^{-4}$ y є = 0.95. ^b1 50 maximum iterations, 20 random values per iteration, $h=1\times {10}^{-4}$ and є = 0.80. ^b2 50 maximum iterations, 20 random values per iteration, $h=1\times {10}^{-4}$ and є = 0.80. ^b3 50 maximum iterations, 20 random values per iteration, $h=1\times {10}^{-4}$ and є = 0.80.

summarizes parameter estimation results across different mesh nodes (7, 15, 21, and 35) using the Luus–Jaakola algorithm. It includes the solution vectors obtained at various operating conditions, their corresponding sum of squares residuals (S), and the CPU time in seconds required for the calculations. Different configurations are highlighted, showing the impact of maximum iterations, random values per iteration, and specific algorithm parameters $(h$ and $ϵ$). The results demonstrate the algorithm’s effectiveness in refining parameter estimates, with solutions indicating reduced residuals and varying computational times based on the node count and settings. This validates the effectiveness of the Luus–Jaakola algorithm in managing complex mesh analyses with various operating conditions.

In this case, the number of iterations influences the precision of the approximation obtained through the Luus–Jaakola algorithm. It is important to highlight that as the number of nodes increases in the numerical solution of PDEs using finite differences, the number of equations to be resolved also rises. Similarly, augmenting the number of iterations in the search process of the Luus–Jaakola algorithm leads to an increase in computational time, as illustrated in Table 2.

Subsequently, the outcomes of parameter estimation utilizing the hybrid approach with varying meshes for spatial coordinate discretization are delineated, along with their respective operational conditions for the hybrid strategy (

Table 3. Parameter estimation results using the hybrid strategy.

Nodes	Solution	Quadratic Sum of Residues (S)	CPU Time(s)
7	[2.96524, 0.97154, 1.98258, 0.50055] a	2.18579 × 10−3	2.0
15	[2.99937, 0.99998, 1.99861, 0.50005] b	2.29720 × 10−6	7.6
21	[2.99969, 0.99996, 1.99943, 0.50001] c	4.20800 × 10−7	41.7
35	[2.99989, 0.99994, 1.99995, 0.49998] d	3.24000 × 10−8	213.0

Operating conditions of the hybrid strategy: ^a Luus–Jaakola (50 maximum iterations, 20 random values per iteration, $h=1\times {10}^{-4}$ and є = 0.95) and Levenberg–Marquardt (λ_ini = 1.0, 22 maximum iterations and $h=1\times {10}^{-4}$); ^b Luus–Jaakola (50 maximum iterations, 20 random values per iteration, $h=1\times {10}^{-4}$ and є = 0.95) and Levenberg–Marquardt (λ_ini = 1.0, 22 maximum iterations and $h=1\times {10}^{-5}$); ^c Luus–Jaakola (50 maximum iterations, 20 random values per iteration, $h=1\times {10}^{-4}$ and є = 0.95) and Levenberg–Marquardt (λ_ini = 1.0, 22 maximum iterations and $h=1\times {10}^{-6}$); and ^d Luus–Jaakola (50 maximum iterations, 20 random values per iteration, $h=1\times {10}^{-4}$ and є = 0.95) and Levenberg–Marquardt (λ_ini = 1.0, 22 maximum iterations and $h=1\times {10}^{-6}$).

). It is noted that parameter estimation with a 35-node mesh achieved the smallest sum of squared residuals (3.24 $\times$ 10⁻⁸), albeit necessitating the longest computation time. This is due to the necessity of solving 99 ordinary differential equations and 6 algebraic equations concurrently at each integration step using the fourth-order Runge–Kutta method.

Figure 4 compares the parameter estimation for a 35-node mesh, including the initial estimation using the Luus–Jaakola algorithm and the subsequent estimation using the Levenberg–Marquardt algorithm.

For the manufactured case studies, the need to implement a hybrid strategy lies in defining an accurate initial approximation close to the actual parameter values when using the Levenberg–Marquardt method. If the initial approximation is inadequate, there is no guarantee of convergence. Therefore, the Luus–Jaakola algorithm should generate an initial approximation that is sufficiently close to ensure convergence of the Levenberg–Marquardt method. In this case, the Luus–Jaakola algorithm provides an appropriate initial approximation that allows the Levenberg–Marquardt method to converge when estimating the four known parameters. If only the Levenberg–Marquardt method is used, it diverges when an initial approximation of 0 or 1 is proposed for each parameter.

3.3. Case Study 3: Controlled Diffusion of an Antimicrobial Peptide in Biopolymer Films

Sebti et al. (2003) [24] conducted an experimental study on the diffusion of nisin in 3% agarose gels. The setup involved a cylindrical agarose gel (70 mm in diameter and height), with its lower face in contact with a large container containing a nisin solution at a 376 μg/mL constant concentration. Initially, the gel was free of nisin, and the solution was slowly stirred, accounting for interfacial resistance. The diffusion coefficient of nisin was estimated from experimental measurements taken at various distances within the gel up to 5.9 days later (Experiment E). Experimental nisin concentrations were determined as a function of distance from the contact face at t = 5.9 days. In this case, we validate the hybrid strategy using the model proposed by Sebti et al. [24] to determine the apparent diffusion coefficient of a nisin solution in contact with an agarose gel, applying Fick’s second law. Figure 5 depicts a one-dimensional schematic diagram of the nisin diffusion experiment within an agarose gel.

Sebti et al. [24] propose that the model for nisin diffusion in agarose gel adheres to Fick’s second law within a semi-infinite cylinder, where the diffusivity coefficient ($D_{eff}$) is constant (Equations (7)–(10)).

$${\displaystyle \frac{\partial C_A}{\partial t}}=D_{eff}{\displaystyle \frac{{\partial }^2{C}_A}{\partial z^2}}$$

(7)

$$\mathrm{B}.\mathrm{C}.1. z=0 C_A=C_A^{\ast }$$

(8)

$$\mathrm{B}.\mathrm{C}.2. z\to \infty { C}_A=C_A^{\ast }=0$$

(9)

$$\mathrm{I}.\mathrm{C}. t=0{ C}_A=C_{init}=0$$

(10)

This problem was analytically resolved using the Laplace transform [24], resulting in the following solution (Equation (11)):

$$W_A={\displaystyle \frac{C_0^{\ast }-C_A}{C_0^{\ast }-C_{init}}}=\mathit{erf}\left({\displaystyle \frac{z}{2\sqrt{D_{eff}t}}}\right)$$

(11)

To validate this case study, experimental data reported by Sebti et al. [24] from Experiment E were obtained using Engauge Digitizer software v. 4.1. The analytical solution was replicated to perform parameter estimation with the hybrid strategy, considering the sum of residuals (Equation (12)):

$$S=\mathit{min}\sum _{i=1}^{data}{\left[W_A^{\mathit{experimental}}-\mathit{erf}\left({\displaystyle \frac{z}{2\sqrt{D_{eff}{t}_c}}}\right)\right]}^2$$

(12)

Additionally, an alternative model proposed by Flores-Martínez et al. [26] accounts for interfacial resistance during the diffusion experiment, assuming that the nisin solution is slowly stirred. This adjustment modifies the original model’s boundary condition, particularly where the agarose gel is in contact with the nisin solution:

$${\displaystyle \frac{\partial C_A}{\partial t}}=D_{eff}{\displaystyle \frac{{\partial }^2{C}_A}{\partial z^2}}$$

(13)

$$\mathrm{B}.\mathrm{C}.1. z=0 {\displaystyle \frac{\partial C_A}{\partial z}}=0$$

(14)

$$\mathrm{B}.\mathrm{C}.2. z=H -D_{eff}{\displaystyle \frac{\partial C_A}{\partial z}}=k_C(C_A-C_{SOL})$$

(15)

$$\mathrm{I}.\mathrm{C}. t=0{ C}_A=C_{init}=0$$

(16)

This model is solved using the method of separation of variables, resulting in the following analytical solution (Equation (17)).

$$C_A=2\left(C_A^0-{\overline{C}}_{Asol}\right){\sum }_{n=1}^{\infty }{\displaystyle \frac{e^{\frac{-{\lambda }_n^2{D}_{\mathit{eff}}t}{H^2}}\mathit{sin}{\lambda }_n}{{\lambda }_n+\mathit{sin}{\lambda }_n\mathit{cos}{\lambda }_n}}\mathit{cos}\left({\lambda }_n\left(1-{\displaystyle \frac{z}{H}}\right)\right)+{\overline{C}}_{Asol}$$

(17)

with the transcendental function (Equation (16))

$${\mathsf{\lambda }}_{\mathrm{n}}\mathrm{s}\mathrm{e}\mathrm{n}{\mathsf{\lambda }}_{\mathrm{n}}-{\mathrm{B}\mathrm{i}}_m\mathrm{c}\mathrm{o}\mathrm{s}{\mathsf{\lambda }}_{\mathrm{n}}=0,\hspace{1em}{Bi}_m={\displaystyle \frac{k_{C}H}{D_{eff}}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{n}=\mathrm{1,2},\mathrm{3,4},5,\dots$$

(18)

This analytical solution was used to validate the numerical solution obtained by applying the hybrid strategy based on the given equations. The model with interfacial resistance was discretized along the spatial coordinate z and solved using the fourth-order Runge–Kutta method on a 51-node mesh. For this case study, knowing approximately the order of magnitude of the parameters to be estimated, an additional scaling factor is added to the Luus–Jaakola algorithm, such that the estimation is achieved with

$$X{p}_i^j=\left(X{p}_i^{j-1}+y_{ki}{r}_i^{j-1}\right)\ast (scaling factor).$$

The models proposed by Sebti et al. [24] and Flores-Martínez et al. [26] were employed to evaluate the hybrid strategy by minimizing the sum of squared residuals, as shown in

Table 4. Results of parameter estimation with the hybrid strategy.

	Analytical Solution (Sebti et al. [24]) a	Analytical Solution (Flores-Martínez et al. [26]) b	Numerical Solution (Flores-Martínez et al. [26]) c
$D_{eff}\left({\displaystyle \frac{{\mathrm{m}}^2}{\mathrm{s}}}\right)$	$8.34183\times {10}^{-11}$	$1.08731\times {10}^{-10}$	$1.08731\times {10}^{-10}$
$K_c\left({\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}\right)$	NA	$1.80632\times {10}^{-8}$	$1.80632\times {10}^{-8}$
$C_A^0\left({\displaystyle \frac{\mu \mathrm{g}}{\mathrm{m}\mathrm{L}}}\right)$	230.8	236.65	236.62
Quadratic Sum of Residues (S2)	2005.89577	1915.81515	1225.28050
MSE	143.27827	136.83394	57.87050
RMSE	11.96989	11.69803	7.60727
CV(RMSE) (%)	14.43087	14.10312	11.10602
Coefficient of Determination (R2)	0.99908	0.99903	0.99944
Average Error	33.28	30.89	27.87

Operating conditions of the hybrid strategy: ^a Luus–Jaakola (200 iterations, 100 random values per iteration and ε = 0.90) and Levenberg–Marquardt (${\lambda }_{init}$ = 1.0, number of iterations = 1000); ^b Luus–Jaakola (100 iterations, 200 random values per iteration and ε = 0.90) and Levenberg–Marquardt (${\lambda }_{init}$ = 1.0, number of iterations = 1000); ^c Luus–Jaakola (50 iterations, 50 random values per iteration, factor1 = 1 × 10⁻⁸, factor2 = 1 × 10⁻⁹, h = 1 × 10⁻⁵ and ε = 0.90) and Levenberg–Marquardt (${\lambda }_{init}$ = 1.0, number of iterations = 1000).

.

Integration of the hybrid strategy yielded an enhanced parameter estimation, leading to a remarkable reduction in the average error rate from 33.28%. This reduction was previously achieved through the utilization of the analytical method based on the Laplace transform elucidated by Sebti et al. [24] to 30.89% via the analytical approach of variable separation explained by Flores-Martínez et al. [26]. Furthermore, the utilization of Equations (13)–(16) was significant in facilitating a more sophisticated numerical estimation, resulting in an average error rate of 27.62%. The sum of squared residuals, MSE, RMSE, CV (RMSE), R^2, and average error statistical metrics show that the numerical solution provided by the hybrid strategy offers a closer fit to the experimental data compared to the analytical solutions, highlighting the significance of accounting for interfacial resistance (see Figure 6).

Flores-Martínez et al. [26] applied the Levenberg–Marquardt technique to perform parameter estimation, resulting in a correlation coefficient of 0.98 between the model and the empirical data. The utilization of a hybrid strategy in parameter estimation through the numerical solution of the model led to a correlation coefficient of 0.99. This establishes the capability of performing parameter estimation utilizing the hybrid strategy even in cases where the analytical solution of a system governed by parabolic-type partial differential equations is challenging to determine.

3.4. Case Study 4. Estimation of Kinetic Parameters in the Conversion of Lactose to Lactic Acid Using Lactobacillus Bulgaricus

In this case, study, we explore the kinetic modeling of lactose conversion to lactic acid using Lactobacillus bulgaricus based on the work of Burgos-Rubio et al. [27]. For this case, the micro-organisms were cultured at 42 °C in a shaker flask for 7 h, with a final inoculum pH of 4.8–5.2, which was placed in a batch-type fermenter using 26.0 g/L of lactose as substrate. The authors proposed different growth kinetics defined by the specific growth rate (μ); these kinetics were adapted to the basic equations that were based on the balance of substrate (S), product (P), and concentration of cellular material (X) using Monod-type kinetics, and the product concentration used was Luedeking–Piret-type kinetics. The biomass, lactic acid, and lactose balances are described in Equations (19)–(22), with the growth rate function defined in Equation (22).

$${\displaystyle \frac{\mathrm{d}\mathrm{X}}{\mathrm{d}\mathrm{t}}}=\mathsf{\mu }\mathrm{X}$$

(19)

$${\displaystyle \frac{\mathrm{d}\mathrm{P}}{\mathrm{d}\mathrm{t}}}=\left(\mathsf{\alpha }\mathsf{\mu }+\mathsf{\beta }\right)\mathrm{X}$$

(20)

$${\displaystyle \frac{\mathrm{d}\mathrm{S}}{\mathrm{d}\mathrm{t}}}=-{\displaystyle \frac{1}{{\mathrm{Y}}_{\mathrm{X}\mathrm{S}}}}\mathsf{\mu }\mathrm{X}-{\displaystyle \frac{1}{{\mathrm{Y}}_{\mathrm{P}\mathrm{S}}}}\left(\mathsf{\alpha }\mathsf{\mu }+\mathsf{\beta }\right)\mathrm{X}$$

(21)

$$\mathsf{\mu }={\mathsf{\mu }}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left({\displaystyle \frac{\mathrm{S}}{\mathrm{S}+{\mathrm{k}}_{\mathrm{s}}}}\right)\left({\displaystyle \frac{K_{PI}}{K_{PI}+P}}\right)$$

(22)

Burgos-Rubio et al. [27] found that the non-competitive product inhibition model (Equation (22)) provided the best fit for the simulation of the kinetics. In this model, seven parameters are estimated, with the following restrictions to ensure physical relevance:

$${\left[{\mathsf{\mu }}_{\mathrm{m}\mathrm{a}\mathrm{x}}\ge 0,{\mathrm{K}}_{\mathrm{s}}\ge 0,{{\mathrm{K}}_{\mathrm{P}\mathrm{I}}\ge 0,\mathsf{\alpha }\ge 0,\mathsf{\beta }\ge 0,\mathrm{Y}}_{\mathrm{X}\mathrm{S}}\ge 0,{\mathrm{Y}}_{\mathrm{P}\mathrm{S}}\ge \mathrm{0,0}\right]}^T$$

(23)

where ${\mathsf{\mu }}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the growth rate (h⁻¹)); ${\mathrm{K}}_s$ is the substrate saturation constant (g/L); ${\mathrm{K}}_{\mathrm{P}\mathrm{I}}$ is the product inhibition constant (g/L)$;{\mathrm{Y}}_{\mathrm{X}\mathrm{S}}$ is the biomass yield coefficient; ${\mathrm{Y}}_{\mathrm{P}\mathrm{S}}$ is the product yield coefficient; α is the growth constant; and $\mathsf{\beta }$ is the non-growth constant of the Luedeking–Piret model.

The system of ordinary differential equations was solved via fourth-order Runge–Kutta method with initial conditions S = 26.0 g/L, P = 0.0 g/L, and X = 0.0009 g/L. Experimental data were digitized using Engauge Digitizer software v. 4.1. from Burgos-Rubio et al. [27]. The initial values for each parameter to estimate the parameters were $\mathrm{X}{\mathrm{p}}_{\mathrm{n}}^{\mathrm{t}}={\left[{\mathsf{\mu }}_{\mathrm{m}\mathrm{a}\mathrm{x}},{\mathrm{k}}_{\mathrm{s}},{\mathrm{Y}}_{\mathrm{X}\mathrm{S}},\mathsf{\alpha },\mathsf{\beta },{\mathrm{Y}}_{\mathrm{P}\mathrm{S}}\right]}^{\mathrm{t}}={\left[1,1,1,1,1,1\right]}^{\mathrm{t}}$. The operating conditions of the hybrid strategy where the lowest sum of squared residuals (7.51369) was obtained were as follows: the estimation with the Luus–Jaakola algorithm was performed over 100 iterations with 50 random values per iteration using a step size for integration (h) of 1 × 10⁻⁴ and a contraction factor (є) of 0.90; the estimation with the Levenberg–Marquardt method used an initial lambda value (λ_init) of 0.1, with an integration step (h) of 1 $\times$ 10⁻⁴ with 10 iterations until reaching convergence, where $\nabla S\le 1\times {10}^{-7}$.

Table 5. Estimated parameters of the ODE system solved by the hybrid approach.

Parameter	Value Obtained in This Work (LJ)	Value Obtained in This Work (LM)	Value Reported by Burgos-Rubio et al. [27]	Units
${\mathsf{\mu }}_{\mathrm{m}\mathrm{a}\mathrm{x}}$	1.78983	1.78976	1.14	h−1
${\mathrm{K}}_{\mathrm{s}}$	1.72357	1.72345	3.36	g/L
${\mathrm{K}}_{\mathrm{P}\mathrm{I}}$	9.31111	9.31089	16.0 a	g/L
$\mathsf{\alpha }$	9.86821	9.86796	9	Dimensionless
$\mathsf{\beta }$	8.39280 × 10−9	0.00000	0	Dimensionless
${\mathrm{Y}}_{\mathrm{X}\mathrm{S}}$	3.67913	3.67913	0.10	g cell/g substrate
${\mathrm{Y}}_{\mathrm{P}\mathrm{S}}$	0.95786	0.95786	0.90	g lactic acid/g substrate
Quadratic Sum of Residues (S2)	7.51369	7.51369	2496.322 a
MSE	0.16697	0.16697	-
RMSE	0.34584	0.34584	-
CV(RMSE)	3.90126	3.90126	-	%
Coefficient of Determination (R2)	0.97250	0.97250	-
Average Error	21.49871	21.49871	-	%

^a Estimated via the Burgos-Rubio et al. [27] approach.

shows the results of the hybrid strategy in micro-organism kinetics.

In this case, the Luus–Jaakola method finds very close approximation to the solution vector found via the Levenberg–Marquardt method, denoting that the quadratic sum of residues is the same in the first five decimal places. Figure 7 illustrates the optimized kinetics of lactose conversion to lactic acid, showing an excellent fit to experimental data, which is critical for scaling up the process from a flask to a bioreactor level.

In all three plots, the close fit between the experimental data and the model predictions underscores the robustness of the hybrid Luus–Jaakola and Levenberg–Marquardt approach to estimating the kinetic parameters. Successfully predicting these key variables is crucial for designing and scaling up the fermentation process from laboratory-scale to industrial-scale bioreactors, ensuring efficient lactic acid production. Additionally, we evaluated four inhibition models to ensure that the non-competitive product inhibition model best satisfies the experimental data. The results are shown in

Table 6. Evaluation of kinetic models in lactic acid production.

Description	QSM	Average Error (%)	MSE	RMSE	CV(RMSE) (%)
Non-Competitive Product Inhibition Model $\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathsf{\mu }={\mathsf{\mu }}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left({\displaystyle \frac{\mathrm{S}}{\mathrm{S}+{\mathrm{k}}_{\mathrm{s}}}}\right)\left({\displaystyle \frac{K_{PI}}{K_{PI}+P}}\right)$	7.51370	21.49871	0.16697	0.34584	3.90126
Competitive Product Inhibition Model 1 a $\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathsf{\mu }={\mathsf{\mu }}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left[{\displaystyle \frac{\mathrm{S}}{\mathrm{S}+{\mathrm{k}}_{\mathrm{s}}(1+\frac{P}{K_{PI}})}}\right]$	12.05869	29.51048	0.26797	0.46181	5.20941
Competitive Product Inhibition Model 2 $\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathsf{\mu }={\mathsf{\mu }}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left({\displaystyle \frac{\mathrm{S}}{\mathrm{S}+{\mathrm{k}}_{\mathrm{s}}}}\right)\left(1-{\displaystyle \frac{P}{K_{PI}}}\right)$	12.49064	28.14847	0.27757	0.48639	5.48671
Uncompetitive Inhibition Model a $\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathsf{\mu }={\mathsf{\mu }}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left[{\displaystyle \frac{\mathrm{S}}{{\mathrm{k}}_{\mathrm{s}}+S(1+\frac{P}{K_{PI}})}}\right]$	8.38326	22.67427	0.18629	0.36403	4.10646

^a Not included in Burgos-Rubio et al. [27].

.

The non-competitive inhibition model exhibits the lowest MAE, MSE, and RMSE values, along with the lowest CV (RMSE) (3.90126%), suggesting that it provides the best fit for the experimental data. Although all three models yield similar R² values (0.97250), the non-competitive model consistently outperforms the others in minimizing residuals and predicting the kinetic behavior of the system. These results confirm that the non-competitive inhibition model remains the most appropriate choice for describing lactic acid inhibition in this system.

3.5. Case Study 5: Parameter Estimation at Heat and Mass Transfer in Drying of Cylindrical Quince Slices

This case is based on the work of Tempezelikos et al. (2015) [25], which investigates the convective drying process of cylindrical quince slices, focusing on the simultaneous heat and mass transfer mechanisms during drying. The study involved developing a numerical model to simulate drying kinetics, temperature distribution, and moisture content within the quince slices. The system geometry is depicted in Figure 8, where drying air flows over a quince slice with thickness L.

The model utilizes coupled partial differential equations to represent heat conduction and moisture diffusion within the slices, considering only the axial direction to capture the cylindrical geometry (L <<< Slice diameter). The boundary conditions are set to include the convective heat and mass transfer at the surface, influenced by drying air temperature and velocity. An understanding of the dynamics of moisture removal and heat distribution within the slices is fundamental to developing effective drying procedures. The primary objective of this case study is to determine the appropriate parameters that govern these mechanisms, including the effective diffusion coefficient and the heat and mass transfer coefficients.

The mechanisms involved in moisture removal from quince are assumed to include the transport of moisture through the product material via liquid diffusion and the evaporation of water from the surface. It is assumed that no phase change, such as evaporation, occurs inside the quince, which is a suitable assumption for diffusion-based models in terms of mass conservation. The driving force for water diffusion inside the quince is the concentration gradient created by the evaporation of water on the surface due to the difference in water vapor partial pressure between the surface and the air. The diffusion of liquid water inside the product is assumed to follow Fick’s law. Hot air supplies heat energy to the quince, part of which is utilized by the evaporation of the water leaving the surface, while the rest is transferred through the material, increasing its temperature.

By integrating these assumptions and understanding the underlying principles, the study aims to develop a robust model for the drying process of cylindrical quince slices. This model will help accurately estimate the parameters necessary for optimizing the drying process, thus ensuring better product quality and energy efficiency. The mathematical model is described as follows in Equations (24)–(35).

Microscopic heat balance

$${\rho }_s{C}_p{\displaystyle \frac{\partial T}{\partial t}}=k_{eff}{\displaystyle \frac{{\partial }^{2}T}{\partial z^2}}$$

(24)

$$\mathrm{B}.\mathrm{C}.1 z=0 +k_{eff}{\displaystyle \frac{\partial T}{\partial z}}=h_{c_1}\left(T-T_{air}\right)+{\lambda }_{\nu }{\rho }_a{k_y}_1(Y_i-Y)$$

(25)

$$\mathrm{B}.\mathrm{C}.2 z=L -k_{eff}{\displaystyle \frac{\partial T}{\partial z}}=h_{c_2}\left(T-T_{air}\right)-{\lambda }_{\nu }{\rho }_a{k_y}_2(Y_i-Y)$$

(26)

Microscopic moisture balance

$${\displaystyle \frac{\partial X}{\partial t}}=D_{eff}{\displaystyle \frac{{\partial }^{2}X}{\partial z^2}}$$

(27)

$$\mathrm{B}.\mathrm{C}.3 z=0 +D_{eff}{\rho }_s{\displaystyle \frac{\partial X}{\partial z}}={k_y}_1{\rho }_a\left(Y_i-Y\right)$$

(28)

$$\mathrm{B}.\mathrm{C}.4 z=L -D_{eff}{\rho }_s{\displaystyle \frac{\partial X}{\partial z}}={k_y}_2{\rho }_a(Y_i-Y)$$

(29)

where $z\in \left[0,L\right]$; L is the axial coordinate, and $t$ is the drying time, ${\rho }_{s}X=C_s$, ${\rho }_{a}Y=C_a$

$$\mathrm{I}.\mathrm{C}. T\left(z,t=0\right)=T_{in} \mathrm{and} X(z,t=0)=X_{in}$$

(30)

${h_c}_1, {h_c}_2$ are the convective heat transfer coefficients at z = 0 and z = L, respectively (W/m² K). ${k_y}_1,{k_y}_2$ are the convective mass transfer coefficients at z = 0 and z = L, respectively (m/s); ${\rho }_s$ is the dry solid density (kg/m³); ${\rho }_a$ is the dry air density (kg/m³); X is the solid moisture on a dry basis (kg H₂O/kg dry solid); and Y is the air humidity on a dry basis (kg H₂O/kg dry air).

The second term on the right-hand side of Equations (25) and (26) (boundary conditions for the heat equation) denotes the superficial moisture vaporization. This value is established by the latent heat of vaporization of water, identified as ${\lambda }_{\nu }$ (J/kg), as described as a Watson equation [25].

$${\lambda }_{\nu }=\mathrm{2,501,050}{\left[{\displaystyle \frac{647.3-\left(T+273.15\right)}{647.3-273.15}}\right]}^{0.3298}$$

(31)

where $p_v^0$ (Pa) is the water vapor pressure [28]

$$p_v^0=9.86923\times {10}^{-5} \mathrm{e}\mathrm{x}\mathrm{p}\left[18.3036-{\displaystyle \frac{3816.44}{T+273.15-46.13}}\right]$$

(32)

$a_{\omega }$ is the water activity of quince [25]

$$a_{\omega }=\mathrm{e}\mathrm{x}\mathrm{p}(-0.1623-3.4766\mathrm{e}\mathrm{x}\mathrm{p}(-7.0063 M\left)\right)$$

(33)

$Y_i$ is the interfacial air humidity [28]

$$Y_i={\displaystyle \frac{18a_{\omega }{p}_v^0}{29(p-a_{\omega }{p}_v^0)}}$$

(34)

$$D_{eff}=D_0 \mathrm{e}\mathrm{x}\mathrm{p}\left[-{\displaystyle \frac{E_a}{R(T+273.15)}}\right]$$

(35)

The effective diffusivity $D_{eff}$ (m²/s) is a function of temperature in an Arrhenius-type model, where $D_0$ is the pre-exponential coefficient (m²/s), and $E_a$ is the activacion energy (J/mol).

The thermophysical properties of quince, as a function of moisture content on a wet basis $M=X/(1+X)$, are as follows: thermal conductivity (W/m K) $k_{eff}=0.148+0.493M$ and specific heat (J/kg K) $c_p=\left(1260+2970M\right)$ [25].

To proceed with the parameter estimation, experimental data on humidity X and temperature T in the solid were digitized, for a velocity and temperature of the drying air of 1 m/s and 40 °C, respectively, using Engauge Digitizer v 4.1.

Table 7. Thermodynamical properties of quince and air.

Property	Value
Slice thickness (L)	0.0092 m
Dry average solid density $\left({\rho }_s\right)$	847.83 kg/m3
Dry air density ${(\rho }_{air})$	1.11 kg/m3
Temperature of drying air $\left(T_{air}\right)$	40 °C
Atmospheric pressure $\left(p\right)$	1 atm
Air humidity ($Y)$	0.00895 kg H2O/kg dry air

shows the constant values used [25]. Due to the geometry of drying, it is expected that the convective coefficients of heat and mass transfer will be greater at z = 0 than at z = L. Tzempelikos et al. calculated them using CFD at various times [25]. Also, the authors calculated effective diffusivity $D_{eff}$ via the slope method, which can be considered as an approximate technique.

Therefore, the indetermined parameters to be estimated are $D_0, E_a,{{ h}_c}_1, {h_c}_2 {k_y}_1,{k_y}_2$. The experimental values of temperature T and dry-basis moisture X were the average values of the corresponding values T = f(z, t) and X = g(z, t); therefore, it is necessary to integrate them from z = 0 to z = L before applying the least squares method [26]. Then, the governing equations were spatially discretized with orthogonal collocation with Legendre polynomials, using 11 interior nodes. This discretization technique integrates variables T and X like a differential quadrature. Then, a trapezoidal-implicit discretization is performed in the time derivative, yielding a set of non-linear algebraic equations that are solved by non-linear relaxation for the Luus–Jaakola approach or Newton–Raphson with LU factorization for the Levenberg–Marquardt approach [26,28,29]. Due to the position of the quince slice in the experimental system, it is expected that the heat and mass transfer coefficients will not be equal at the inlet and outlet of the air flow.

Table 8. Results of parameter estimation in the drying of quince ^a.

Parameter	Value	Unity
$D_0$	4.62926 × 10−4	m2/s
${\displaystyle \frac{E_a}{R}}$	3966.07274	K
${{ h}_c}_1$	241.83632	W/m2 K
${h_c}_2$	10.57547	W/m2 K
${k_y}_1$	0.64438	m/s
${k_y}_2$	2.60012 × 10−3	m/s
Quadratic Sum of Residues (S2)	4.04737
MSE	0.08799
RMSE	0.23003
CV(RMSE)	1.30688	%
Coefficient of Determination (R2)	0.99983
Average Error	2.11108	%

^a Tzempelikos et al. [25].

shows the results of applying the Luus–Jaakola method (200 iterations, 50 random values per iteration and ε = 0.90) and the Levenberg–Marquardt method (${\lambda }_{init}$ = 1.0, number of iterations = 20 and residual = 10⁻⁵).

Figure 9 shows an excellent agreement of the predictions of the governing equations with the estimated parameters, comparing the experimental values of temperature and moisture of the quince slices, in concordance with the statistical metrics obtained. A mesh independence analysis was performed with 15 and 19 interior orthogonal collocation points, which validated the mesh size used. Therefore, parameter estimation is the activity before scaling and designing processes from the governing equations. In this case, the Levenberg–Marquardt method requires a very close approximation to the solution vector, since otherwise, there is no convergence, or oscillating convergence, due to the high non-linearity of the governing equations. Thus, the Luus–Jaakola method is decisive in finding a close approximation that allows LM to reach convergence. This case study emphasizes the hybrid strategy’s capacity to handle the complexities of simultaneous heat and mass transfer in drying processes, providing accurate parameter estimates essential for optimizing industrial drying operations.

4. Conclusions

This study successfully demonstrates the effectiveness of a hybrid optimization strategy that combines the strengths of the Luus–Jaakola (LJ) and Levenberg–Marquardt (LM) algorithms for complex parameter estimation tasks. The LJ algorithm, known for its robustness in global tasks, was particularly effective in generating accurate initial approximations, which are crucial for the LM algorithm to converge on precise solutions in non-linear and multi-parameter systems.

The hybrid approach was tested across various case studies, each involving different types of differential equations, including parabolic partial differential equations (PDEs), ordinary differential equations (ODEs), and algebraic equations. The results consistently showed that the LJ algorithm could navigate the complex landscape of these equations to find multiple local minima. This capability, combined with the LM algorithm’s efficiency in refining these initial estimates, ensured that the method could reliably achieve convergence to an optimal solution.

One of the key strengths of this hybrid strategy is its ability to address the challenges posed by parameter estimation in systems with significant variations in parameter magnitudes, such as those found in heat and mass transfer processes. By ensuring that the initial approximation is close to the solution, the hybrid method reduces the risk of divergence and improves the accuracy of the final parameter estimates.

Overall, this hybrid LJ-LM approach has proven to be a versatile and powerful tool for solving complex parameter estimation problems across a range of scientific and engineering disciplines. Its ability to overcome the limitations of traditional methods by combining global and local search strategies makes it particularly valuable for applications involving complex mathematical models and systems governed by differential equations.