Identification of Crude Distillation Unit: A Comparison between Neural Network and Koopman Operator

Abstract

In this paper, we aimed to identify the dynamics of a crude distillation unit (CDU) using closed-loop data with NARX−NN and the Koopman operator in both linear (KL) and bilinear (KB) forms. A comparative analysis was conducted to assess the performance of each method under different experimental conditions, such as the gain, a delay and time constant mismatch, tight constraints, nonlinearities, and poor tuning. Although NARX−NN showed good training performance with the lowest Mean Squared Error (MSE), the KB demonstrated better generalization and robustness, outperforming the other methods. The KL observed a significant decline in performance in the presence of nonlinearities in inputs, yet it remained competitive with the KB under other circumstances. The use of the bilinear form proved to be crucial, as it offered a more accurate representation of CDU dynamics, resulting in enhanced performance.

1. Introduction

Crude distillation units (CDUs) play a significant role in the refining sector, responsible for fractionating crude oil into intermediate products that are subsequently processed in downstream units to meet the market specifications. The quality of these intermediate products, primarily influenced by the operating conditions of the CDUs, is crucial for ensuring the quality of the final refinery products [1]. Currently, the oil refining industry is confronting several challenges, including a significant increase in crude oil prices, fluctuations in product demand driven by market dynamics, and specific regulatory constraints imposed on industrial activities [2]. Improving these units may lead to significant improvements in the efficiency and reliability of oil refining processes, ultimately reducing operational costs and increasing overall savings. Developing suitable process control for the CDU is necessary to achieve these objectives. Nevertheless, the existing body of literature primarily focuses on the design of control algorithms, rather than the identification of a CDU to a sufficient accuracy.

In recent years, significant advancements in artificial intelligence have facilitated the development of dynamic models through various data-driven techniques, such as polynomial regressions [3,4,5,6], support vector regression (SVR) [7], and Artificial Neural Networks (ANNs). For example, Liau et al. [8] aimed at optimizing product outputs by using an ANN to predict the yield of kerosene, diesel, and atmospheric gas oil. Motlaghi et al. [9] employed an ANN to predict products flow rate that were optimized based on market values. Gueddar et al. [10] developed an ANN model to optimize energy efficiency by considering crude oil properties, such as boiling point and crude flow rate. Building on this approach, Durrani et al. [11] developed a multi-output ANN model to address variations in crude composition and predict optimum cut-point temperatures, using a hybrid Taguchi and genetic algorithm for more energy-efficient operations. Ochoa-Estopier et al. [12] developed an ANN model for a CDU and employed a Simulated Annealing (SA) optimizer to enhance revenue while reducing energy usage. This work was further extended by the same authors in [13,14], who incorporated a heat exchanger network model to enhance operational optimization, aiming to boost net profit while adhering to practical constraints. Shi et al. [15] modeled a CDU process using a wavelet neural network which was combined with the line-up competition algorithm (LCA) for the economic optimization of the CDU operation. More recently, a Long Short-Term Memory network (LSTM) was developed by [16] to predict and analyze energy efficiency in the CDU under different operating conditions. A hybrid ANN-SVM model was developed by [17] to simulate the performance of the CDU accurately and efficiently within an optimization framework. Li et al. [18] developed a hybrid Fuzzy Logic–ANN model to construct a knowledge-based strategy to adapt to different feedstock properties. Bootstrap ANN models were used by Osuolale and Zhang in [19] and Muhsin et al. in [20] to develop a model of a CDU process, with the former authors focusing on energy efficiency and the latter on maximizing the production rate. A comparison between different data-driven models for predicting CDU product properties was investigated by [21], including PCA-ResNet, SOM-ResNet, Feedforward Neural Networks (FNNs), Partial Least Squares (PLS), and LASSO. The study concluded that incorporating prior knowledge and employing appropriate dimensionality reduction techniques, such as PCA-ResNet, greatly improved model accuracy. Using machine learning can offer more accurate and sometimes computationally efficient solutions compared to the complex and resource-intensive nature of ANNs [22]. Fadzil et al. [23] explored five machine learning models for optimizing product yields based on varying feed properties and operating conditions. These models included decision tree regression, support vector regression, ANN, random forest regression, and extreme gradient boosting (XGBoost), with XGBoost demonstrating superior performance.

While ANNs have proven to be suitable for modeling a CDU, they rely on static mapping of outputs from inputs using data; this limits our understanding of the physical or chemical mechanisms governing the CDU. Linear predictors can be utilized for nonlinear systems to effectively capture and model nonlinear behavior while also identifying the dynamic response of a linear system. From a control engineering perspective, this approach provides valuable insights into the process by analyzing its time-domain characteristics [24]. In this context, Bernard Koopman [25] introduced his operator, which describes the evolution of measurements in Hamiltonian systems over time. Mezic and his collaborators [26,27] expanded these concepts to include nonconservative systems, extending its utility to a wide range of applications, solar panels [28], power systems [29], robotics [30], autonomous driving [31], biology [32], and traffic flow [33], to list a few.

The quality of a model is usually evaluated by its ability to generalize to new, unseen data. In regression tasks, if the model is properly selected in terms of structure and hyperparameters, and overfitting is avoided, the problem can be considered conceptually solved. However, when dealing with real data from advanced process control units, inputs to the model may fall outside the training domain due to factors like model-plant mismatches, poor tuning, or stringent constraints in the closed-loop system. In such situations, it is essential for the regression model to make acceptable predictions or, at the very least, not fail. This subject is the topic of this paper. In this work, a comparative analysis of modeling a CDU under different experimental conditions is conducted. This includes the Koopman operator in both linear (KL) and bilinear (KB) forms, as well as a NARX−NN model. The performance of these models is tested in real-world settings, such as gain and delay mismatches, nonlinearities, and disturbances. Bayesian optimization is used for hyperparameter tuning to ensure a fair comparison. The remainder of the paper is structured as follows: Section 2 provides preliminaries of the methodologies used. Section 3 covers the process description and data generation. The results and discussion are presented in Section 4. Finally, the conclusion and suggestions for future work are outlined in Section 5.

2. Modeling and Methodology

2.1. Nonlinear Autoregressive Exogenous Model

The Nonlinear AutoRegressive Exogenous model (NARX) was proposed as a novel framework for representing a diverse range of nonlinear systems. The fundamental principle behind this model is the incorporation of past inputs and outputs, which helps to mitigate the effects of both nonlinearities and delays in the output [34]. This model is defined as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y\left(k\right)=F[& y(k-1),y(k-2),\dots ,y(k-n_y),\hfill \\ & u(k-d),u(k-d-1),\dots ,u(k-d-n_u)]\hfill \end{array}\end{array}$$

(1)

where $y\left(k\right)$ and $u\left(k\right)$ represent the output and the input of the system, respectively, with $n_y$ and $n_u$ the maximum lags for the output and input, whereas d is the time delay. $F[.]$ denotes a nonlinear mapping function. In practical applications, there are several methods that can be used to approximate the unknown function $F[.]$; these include polynomial series, wavelet expansions, Fuzzy Logic, Neural Networks, and many others.

While the formulation presented in Equation (1) is designed for a Single-Input Single-Output (SISO) system, it is straightforward to extend these models to accommodate Multi-Input Multi-Output (MIMO) systems. For example, consider a generic MIMO scenario with r inputs $\{u_1\left(k\right),u_2\left(k\right)..,u_r\left(k\right)\}$ and s outputs $\{y_1\left(k\right),y_2\left(k\right)..,y_s\left(k\right)\}$, as given by Equation (2).

$$\begin{array}{c}\hfill \begin{array}{cc}& u_i^{[k-1]}=[u_i(k-1),u_i(k-2),\dots ,u_i(k-n_u)]\hfill \\ & y_j^{[k-1]}=[y_j(k-1),y_j(k-2),\dots ,y_j(k-n_y)]\hfill \\ & i=1,2,3,\dots ,rj=1,2,\dots ,s\hfill \end{array}\end{array}$$

(2)

The formulation of the NARX model for the MIMO system (2) can be expressed as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y_1\left(k\right)& =F_1[y_1^{(k-1)},\dots ,y_s^{(k-1)},u_1^{(k-1)},\dots ,u_r^{(k-1)}]\hfill \\ \hfill y_2\left(k\right)& =F_2[y_1^{(k-1)},\dots ,y_s^{(k-1)},u_1^{(k-1)},\dots ,u_r^{(k-1)}]\hfill \\ \hfill ⋮\\ \hfill y_s\left(k\right)& =F_s[y_1^{(k-1)},\dots ,y_s^{(k-1)},u_1^{(k-1)},\dots ,u_r^{(k-1)}]\hfill \end{array}\end{array}$$

(3)

where $F_1[.]$, $F_2[.]$,..., $F_s[.]$ are nonlinear functions. ANNs are widely used in signal processing, pattern recognition, data fitting, analysis, and control due to their ability to learn and represent mathematical descriptions of systems [35]. This can be either Single-Layer Networks (SLNs) or Multi-Layer Networks (MLNs). SLNs, in particular, have been extensively applied to reconstruct the mapping function $F[.]$, resulting in more accurate NARX models in different applications [36,37,38,39]. Therefore, this approach was adopted in our work. Figure 1 depicts the conventional configuration of a NARX−NN model. Combining the past inputs and outputs into a single vector $\mathbf{x}\left(\mathbf{k}\right)$:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{x}\left(\mathbf{k}\right)=[& u_1(k-1),\dots ,u_1(k-n_u),\hfill \\ & u_r(k-1),\dots ,u_r(k-n_u),\hfill \\ & y_1(k-1),\dots ,y_1(k-n_y),\hfill \\ & y_s(k-1),\dots ,y_s(k-n_y){]}^T\hfill \end{array}\end{array}$$

(4)

The NARX−NN model can be mathematically expressed in a matrix form as follows:

$$\mathbf{y}(\mathbf{k}+\mathbf{1})=\varphi \left(\mathbf{W}\mathbf{x}\right(\mathbf{k})+\mathbf{b})$$

(5)

where $\mathbf{y}(\mathbf{k}+\mathbf{1})$ is the output vector $\{y_1(k+1),\dots ,y_s(k+1)\}$, $\mathbf{W}$ is the weight matrix, $\mathbf{b}$ is the bias vector, and $\varphi$ is typically a nonlinear activation function applied element-wise (e.g., sigmoid, tanh, ReLU). A pseudo-code for the identification of the NARX−NN model is presented in Algorithm 1.

Algorithm 1 NARX−NN model identification algorithm

1: Input Data $D=Y,U$  2: Initialize $\mathbf{W},\mathbf{b}$  3: for all epochs do  4:     for all batches of D do  5:     Construct input vector from past inputs/outputs Equation (4)  6:      Forward pass to obtain predicted output $\widehat{\mathbf{y}}(\mathbf{k}+\mathbf{1})$ Equation (5)  7:         Compute loss $\mathcal{L}={\displaystyle \frac{1}{N}}{\sum }_{k=1}^N{(\mathbf{y}(\mathbf{k}+\mathbf{1})-\widehat{\mathbf{y}}(\mathbf{k}+\mathbf{1}))}^2$  8:         Update $\mathbf{W}$ and $\mathbf{b}$ using backpropagation  9:     end for 10: end for

2.2. Koopman Operator Theory

In this section, a brief introduction to the Koopman operator is provided; for more details, readers should refer to [40].

First, consider a discrete nonlinear autonomous system in the following form:

$$x_{k+1}=f\left(x_k\right)$$

(6)

where $x_k\in \mathbb{X}\subset {\mathbb{R}}^n$ is the vector of state variables, $\mathbb{X}$ is the state space, and $f:\mathbb{X}\to \mathbb{X}$ is a map function that represents the evolution of the state. The infinite-dimensional Koopman operator $\mathcal{K}$ acting on a set of observation functions g is defined as

$$\mathcal{K}g\left(x_k\right)=g\left(f\left(x_k\right)\right)=g\left(x_{k+1}\right)$$

(7)

where $\mathcal{K}:\mathbb{G}\to \mathbb{G}$ and $g:\mathbb{X}\to \mathbb{G}$. $\mathbb{G}$ denotes the lifted space. The Koopman operator advances the value of the observable function forward in time. One might obtain a finite approximation of the Koopman operator by restricting it to an invariant subspace ${\mathbb{G}}_{\mathbb{N}}\subset \mathbb{G}$, which is spanned by a set of Koopman eigenfunctions that satisfy the following property

$$\lambda \varphi \left(x\right)=\phi \left(f\left(x\right)\right)$$

(8)

where the eigenvalue $\lambda$ corresponds to the eigenfunction $\phi \left(x\right)$. Using the Koopman Mode Decomposition (KMD) [27], the observation function $g\left(x_k\right)$ can be represented as

$$g\left(x\left(k\right)\right)={\mathcal{K}}^{k}g\left(x_0\right)=\sum {\lambda }_j^k{\phi }_j\left(x_0\right)v_j$$

(9)

where $v_j$ is the jth Koopman mode associated with the eigenfunction ${\phi }_j$ and the eigenvalue ${\lambda }_j$.

For a discrete-time control system, where $u\in \mathbb{U}\subset {\mathbb{R}}^m$,

$$x_{k+1}=f(x_k,u_k)$$

(10)

We define an extended state $\tilde{x}:={[x^T,u^T]}^T$, and let $g(x,u):\mathbb{X}\times \mathbb{U}\to \mathbb{G}$ be a set of scalar observable function of the extended state space. The Koopman operator advances measurement functions according to

$$\mathcal{K}g(x_k,u_k)=g(f(x_k,u_k),u_k)=g(x_{k+1},u_{k+1})$$

(11)

The Koopman eigenpairs $(\phi ,\lambda )$ associated with Equation (11) satisfy

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{K}\phi (x_k,u_k)& =\phi (f(x_k,u_k),u_k)\hfill \\ & =\phi (x_{k+1},u_{k+1})=\lambda \phi (x_k,u_k)\hfill \end{array}\end{array}$$

(12)

Assuming that the observation function can be separated into two parts:

$$g\left(x,u\right)=\left[g_x\left(x,u\right);g_u\left(x,u\right)\right]$$

(13)

Then, the expression for the Koopman operator would have the following form [41,42]:

$$\left[\begin{array}{c}{g}_x\left(x_{k+1},u_{k+1}\right)\\ u_{k+1}\end{array}\right]=\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]\left[\begin{array}{c}{g}_x\left(x,u\right)\\ g_u\left(x,u\right)\end{array}\right]$$

(14)

where $A,B,C,D$ are submatrices of the Koopman operator matrix $\mathcal{K}$. We assume that $g_x\left(x,u\right)$ depends only on the system’s state, $g_x\left(x,u\right)=g_x\left(x\right)$. An embedding ANN, $\vartheta \left(x_k\right):{\mathbb{R}}^n\to {\mathbb{R}}^d$, where d is a tunable hyperparameter, shall be used to learn the observation function $g_x\left(x\right)$. Furthermore, the control inputs in the lifted space are assumed to remain unchanged compared to those in the nonlinear space; hence, $g_u\left(x,u\right)=u$. Denoting $z_k=g_x\left(x_k,u_k\right)$ as the state within the lifted dynamics, a linear relationship of the state variable z can be obtained as follows

$$z_{k+1}=A{z}_k+B{u}_k$$

(15)

By employing the Koopman Canonical Transformation (KCT) [43], the Koopman bilinear form (KBF) [44] is given by,

$$z_{k+1}=A{z}_k+\sum _{i=1}^m{B}_i{u}_i$$

(16)

A schematic representation of the Koopman dynamics based on an ANN is presented in Figure 2.

Furthermore, the original states are concatenated with the output of $\vartheta$, allowing for the direct retrieval of the original states:

$$\begin{array}{cc}\hfill & z_k=\left[\begin{array}{c}{x}_k\\ \vartheta \left(x_k\right)\end{array}\right]\hfill \\ \hfill & x_k=C_x{z}_k,C_x=\left[\begin{array}{cc}{I}_n& 0\end{array}\right]\hfill \end{array}$$

(17)

To learn the Koopman dynamics jointly with the embedding network, we formulate the following loss function

$$\begin{array}{cc}\hfill \mathcal{L}=& {\mathcal{L}}_{pred}+\alpha {\mathcal{L}}_{lin}+\beta {\mathcal{L}}_{L1}\hfill \end{array}$$

(18)

where

$$\begin{array}{cc}\hfill {\mathcal{L}}_{pred}=& {\parallel x_{k+1}-C_x{\widehat{z}}_{k+1}\parallel }^2\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\mathcal{L}}_{lin}=& {\parallel z_{k+1}-{\widehat{z}}_{k+1}\parallel }^2\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\mathcal{L}}_{L1}=& {\parallel A\parallel }_1+{\parallel B\parallel }_1\hfill \end{array}$$

(21)

where ${\parallel .\parallel }^2$ represents the Mean Squared Error (MSE), averaged over the data points, then over the number of outputs. The term $\widehat{z}$ denotes the estimated lifted state obtained from Equation (15) if a linear model is used, or from Equation (16) if a bilinear model is employed. The loss term ${\mathcal{L}}_{pred}$ represents the one-step prediction error in the nonlinear state-space, which ensures the consistency of the Koopman model with the underlying nonlinear system as it progresses over time. ${\mathcal{L}}_{lin}$ represents the one-step prediction error in the lifted space. Additionally, ${\mathcal{L}}_{L1}$ encourages sparsity in the Koopman dynamical system, facilitating better generalization and reducing overfitting. The coefficients $\alpha$ and $\beta$ determine the weighting of ${\mathcal{L}}_{lin}$ and ${\mathcal{L}}_{L1}$ with respect to ${\mathcal{L}}_{pred}$. The algorithm for training Koopman dynamics is summarized as Algorithm 2, where ${\vartheta }_{\theta }$ are the weights of the embedding network $\vartheta$.

Algorithm 2 Training procedure of the Koopman dynamics

Input Dataset $D=X,U$ Initialize $A,B,{\vartheta }_{\theta }$ for all epochs do     for all batches of D do         lift the state $z_k\leftarrow [x_k,\vartheta \left(x_k\right)]$         lift the state $z_{k+1}\leftarrow [x_{k+1},\vartheta \left(x_{k+1}\right)]$        Forward $z_k$ one step in time to obtain the estimated lifted state $\widehat{z}$         Compute prediction loss Equation (19)         Compute linear loss Equation (20)         Compute the total loss Equation (18)         Update $A,B,{\vartheta }_{\theta }$ by back propagation     end for end for

3. Application to CDU

3.1. Process Description

The crude distillation unit (CDU) is the initial processing unit in a petroleum refinery, responsible for separating crude oil into various distillate streams that serve as essential raw materials for downstream refining processes. This process can be divided into several units based on the level of output, including the pre-flash unit (PFU), atmospheric distillation unit (ADU), vacuum distillation unit (VDU), splitting unit (SPU), stabilizer unit (SBU), and heat exchanger network (HEN). Figure 3 illustrates a basic representation of the ADU used in this work. In the CDU, crude oil is first heated, and water is injected to dissolve salts. The mixture then undergoes electrostatic precipitation in a desalter drum to separate the salts. The crude oil is subsequently routed to the distillation tower’s flash zone, where it is heated in a fired heater to vaporize distillate products. Overflash is added to ensure effective reflux streams within the tower. The heated crude oil enters the fractionation tower in the flash zone, where the distillate vapors ascend the tower, countering a colder liquid reflux stream. Distillate products are then removed from selected trays, stream stripped, and sent to storage [45,46].

3.2. Data Collection

The data were generated using a simulated CDU unit under a Model Predictive Controller (MPC) within the Aspen DMC Plus software. The data included three control variables (CVs), three manipulated variables (MVs), and one feedforward (FF) variable, as listed in

Table 1. List of inputs and outputs.

Inputs
$u_1$	Column reflux
$u_2$	Middle draw
$u_3$	Column bottom tray temperature
d	Feed to column
Outputs
$y_1$	Overhead composition
$y_2$	Middle-draw composition
$y_3$	Bottoms’ composition

. The sampling period was set to 1 min, with a time to steady state of 30 min. A prediction horizon of 30 and a control horizon of 10 were used for the MPC. Various experimental conditions were simulated to depict real-life situations, including gain, delay, and time-constant mismatches, tight constraints, nonlinearities, and poor PID tuning. Each run lasted for 1440 samples, equivalent to 1 day.

Table 2. Description of each experimental condition.

Case #	Description
1	No mismatch with aggressive tuning
2	Gain mismatch, $u_1-y_1,u_1-y_2$
3	Gain mismatch, $u_2-y_2$
4	Delay mismatch, $u_1-y_1,y_2,y_3$
5	Poor PID tuning for all three MVs
6	Sqrt-type nonlinearity in $u_3-y_1,y_2,y_3$ channels
7	Tight constraints on $u_1,u_2,u_3$
8	Unmeasured disturbance affecting $y_1,y_2,y_3$ equally

describe the settings for each case. To validate the efficiency of the proposed methods in capturing the dynamics of the CDU, training and validation were performed on the nominal dataset (Case 1), and each model was tested on the remaining datasets.

4. Results and Discussion

To ensure a fair comparison among the models, the hyperparameters of each model were fine-tuned using Bayesian optimization, implemented with the KerasTuner API. The objective was to minimize the validation MSE.

Table 3. Hyperparameter search space.

Model	Hyperparameter	Range
KL / KB	Network width	[16, 256]
	Network depth	[1, 32]
	No. of lifted dynamics	[4, 64]
NARX	Network width	[4, 256]
	Network depth	[1, 32]
	Delay	[1, 4]

outlines the specific hyperparameters and their respective search ranges. The results are summarized in

Table 4. Hyperparameters tuning results.

Hyperparameter	NARX	KL	KB
Network width	12	16	22
Network depth	1	8	10
Delay	2	/	/
No. of lifted dynamics	/	42	64

. In what follows, the Koopman linear and bilinear models are referred to as KL and KB, respectively.

The results based on the training dataset are summarized in Figure 4, and those for the test datasets are provided within Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11. The plots in Figure A1, Figure A2, Figure A3, Figure A4, Figure A5, Figure A6, Figure A7, Figure A8, Figure A9, Figure A10, Figure A11, Figure A12, Figure A13, Figure A14, Figure A15, Figure A16, Figure A17, Figure A18, Figure A19, Figure A20, Figure A21, Figure A22, Figure A23 and Figure A24 in Appendix A depict the predicted outputs ($y_1$, $y_2$, $y_3$) vs. the actual values for each case. As can be noticed from Figure 4, the NARX−NN model had the lowest average MSE (0.0081) on the training dataset, suggesting that NARX−NN fitted the training data well and indicating good initial performance in the identification of the CDU. In comparison, the KB model performed slightly better than the KL model, where the latter showed the highest average MSE of 0.0127. However, across all test scenarios, KB consistently achieved superior performance with the lowest average MSE, indicating better generalization ability and robustness in capturing the dynamics of the CDU across different experimental conditions. The KL model performed competitively with KB and often outperformed NARX−NN. While benefiting from exogenous inputs and past outputs is powerful for NARX−NN, it also adds more complexity to the model. This may be sensitive in terms of model parameters and can lead to overfitting, causing poorer generalization in capturing the dynamics under different experimental conditions.

In cases where nonlinearities were present in the inputs, both the KL and KB models experienced a noticeable degradation in performance. This degradation was primarily due to the assumption that the control inputs remained unchanged within the lifted dynamics. However, the KL model exhibited the poorest performance; conversely, the inclusion of a bilinear term in the KB model allowed for a more flexible representation of the system dynamics. This flexibility enabled the KB model to capture certain nonlinear behaviors that the linear KL model could not, leading to its superior performance. Based on these observations, one may conclude that under different experimental conditions, the NARX−NN model dramatically failed to provide predictions, while the KB model continued to make good predictions.

5. Conclusions

This article presented a comparison of three different models in identifying the CDU process, namely NARX−NN, KL, and KB, conducted under different experimental conditions, such as gain, delay mismatches, nonlinearities, and disturbances. During the training process, the NARX−NN model had the lowest average Mean Squared Error (MSE), indicating a high level of initial competency. Nevertheless, KB demonstrated higher performance and robustness across different datasets used for testing, surpassing the performance of both NARX−NN and KL. Although KL showed comparable performance to KB and sometimes outperformed NARX−NN, its performance significantly deteriorated when faced with nonlinearities.

The use of a bilinear term in KB was the main differentiating factor, providing a more flexible representation of system dynamics compared to the linear assumption in the KL model. This flexibility allowed KB to capture nonlinear phenomena that other models could not, resulting in superior performance, as evidenced by lower MSE values. The KB stands out as a potential option because of its improved capacity to generalize.