Generalized bilinear Koopman realization from input-output data for multi-step prediction with metaheuristic optimization of lifting function and its application to real-world industrial system

1 Generaliz ed bilinear K oopman realization from input-output data f or multi-step prediction with metaheur istic optimization of lifting function and its application to real-world industr ial system Shuichi Y ahagi, Membership , Ansei Y onezaw a, Membership , Heisei Y onezawa, Membership , Hiroki Seto , Membership , and Itsuro Kajiwara Abstract —This paper introduces an input-output bilin- ear K oopman realization with an optimization algorithm of lifting functions. For nonlinear systems with inputs, K oopman-based modeling is effective because the K oop- man operator enables a high-dimensional linear repre- sentation of nonlinear dynamics. Howe ver , traditional ap- proac hes face significant challenges in industrial appli- cations. Measuring all system states is often impractical due to constraints on sensor installation. Moreover , the predictive perf ormance of a K oopman model strongly de- pends on the choice of lifting functions, and their design typically requires substantial manual effort. In addition, although a linear time-in variant (L TI) Koopman model is the most commonly used model structure in the K oopman framew ork, such model exhibit limited predictive accuracy . T o address these limitations, we propose an input-output bilinear Koopman modeling in which the design parame- ters of radial basis function (RBF)-based lifting functions are optimized using a global metaheuristic algorithm to impro ve long-term prediction performance. Consideration of the long-term prediction performance enhances the reli- ability of the resulting model. The proposed methodology is validated in simulations and experimental tests, with the airpath control system of a diesel engine as the plant to be modeled. This plant represents a challenging industrial application because it exhibits strong nonlinearities and coupled multi-input multi-output (MIMO) dynamics. These results demonstrate that the proposed input-output bilinear K oopman model significantly outperforms traditional linear K oopman models in predictive accuracy . Index T erms —Bilinear system, Data-driven modeling, K oopman operator , Nonlinear system, Diesel engine, Indus- trial application Manuscript received Month xx, 2xxx; revised Month xx, xxxx; ac- cepted Month x, xxxx. This work was suppor ted in par t by the ISUZU Advanced Engineer ing Center , Ltd. Shuichi Y ahagi (Corresponding author) is with Depar tment of Mechan- ical Engineer ing, T okyo City University , 1-28-1 T amazutsumi, Setagay a- ku, T okyo 158-8557, Japan (e-mail: yahagisi@tcu.ac.jp). Ansei Y onezaw a is with Depar tment of Mechanical Engineer ing, K yushu University , 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan. Heisei Y onezawa and Itsuro Kajiwara are with Division of Mechanical and Aerospace Engineering, Hokkaido University , N13, W8, Kita-ku, Sapporo , Hokkaido 060-8628, Japan. Hiroki Seto are with 6th Research Depar tment, ISUZU Advanced Engineering Center Ltd., 8 Tsutidana, Fujisawa-shi, Kanagawa 252- 0881, Japan. I . I N T RO D U C T I O N A. Motiv ation S Ystem identification is essential for enabling system pre- diction, control, and fault detection. Traditional system identification methods, such as Auto-Re gressi ve with eX- ogenous (ARX) models and Numerical Algorithms for Sub- space State Space System Identification (N4SID) [1] hav e prov en highly effecti ve for linear systems. Howe ver , identi- fying models for highly nonlinear industrial systems remains challenging, which in turn complicates the control of complex industrial processes. T o address this issue, data-driv en ap- proaches for nonlinear systems hav e been proposed, including the Koopman operator framework [2], [3], [4], deep learning methods such as recurrent neural networks (NNs) and long short-term memory (LSTM) networks [5], nonlinear ARX (N ARX) models [6], and sparse identification of nonlinear dynamics (SINDy) [7], [8]. The Koopman operator is formulated as a linear operator in an infinite-dimensional space by mapping the original state space to a higher-dimensional space through lifting functions [9], [10]. Its ability to represent a nonlinear sys- tem as a linear model is a distincti ve feature not found in other data-dri ven approaches. Although Koopman operator theory inherently deals with infinite-dimensional spaces, in practical computations, a finite-dimensional approximation of the K oopman operator can be identified from data using dynamic mode decomposition (DMD) [11], [12] and extended DMD (EDMD) [13], [14]. The resulting K oopman model can accurately capture the original nonlinear behavior by le ver- aging a high-dimensional linear structure. This enables the application of con ventional linear control theory [2], [15], [16]. Despite these achievements, sev eral challenges for application to industrial systems remain, including the selection of lifting functions that significantly impact prediction performance, the inherent limitations of linear model predictions, and the unav ailability of measurements for the full system states. B. Related work It has been noted that L TI Koopman models may not adequately capture the control-affine dynamics of nonlinear systems, which makes accurate prediction challenging [17]. 2 T o address this limitation, bilinear K oopman realizations have been inv estigated. These models of fer improved predictiv e accuracy compared to L TI K oopman models while remaining computationally more ef ficient than fully nonlinear K oopman approaches [18]. Theoretical foundations of bilinear Koopman realizations have been explored in [18], [19]. In addition, methods for computing bilinear lifting functions using deep NNs have been proposed [20], [21], [22]. The choice of the observation function that maps state variables to observable variables–commonly referred to as the lifting function–has a significant impact on the predicti ve performance of the K oop- man model. According to [23], previous approaches for de- signing lifting functions include mechanics-inspired selections, empirical selections such as monomials and polynomials, and RBF functions with randomly assigned centers. Another pro- posed approach employs highly contributory basis functions deriv ed from SINDy modeling [24], [25]. Howe ver , these methods often fail to ensure reproducibility or improve the predictiv e accuracy of model identification. T o address these limitations, recent research has explored neural network-based approaches for lifting function design [26], [27]. Nevertheless, NN-based methods present challenges such as hyperparameter tuning, overfitting, model complexity , computational cost, and issues related to vanishing or exploding gradients. Further- more, the orthogonality of NN-based lifting functions is not guaranteed. T o date, meta-heuristic approaches such as particle swarm optimization (PSO) ha ve not been inv estigated for lifting function design. Previously , L TI and bilinear Koopman realizations typically assumed that all system states and inputs were measurable. Howe ver , in practical industrial applications, it is often infea- sible to deploy sensors capable of capture every state variable. T o address this limitation, a Koopman realization framework that relies solely on input-output data is required. This per- spectiv e is particularly important for industrial implementa- tion. Moreover , prior research on K oopman operators using input-output data [2], [25], [28] has been limited, with most studies focusing exclusi vely on L TI Koopman formulations. Additionally , the selection of the arguments (i.e., embedded states) of the lifting function for controlled systems has not been thoroughly in vestigated. Although modeling of diesel engine airpath systems has been widely in vestigated [29], [30], [31], achie ving accurate and robust models remains highly challenging. These systems exhibit nonlinear and multiv ariable dynamics. The traditional approaches rely on physical modeling [30], often represented as linear parameter -varying (LPV) systems [32]. Howe ver , physical modelings suf fer from limited accuracy and dif ficul- ties in parameter identification. In data-dri ven approaches, L TI modeling has been performed using ARX identification [33]; howe ver , the applicability of L TI models is restricted due to the wide operating range of engines. For nonlinear system identification, studies such as [34] have adopted NARX-based methods, while the literature [29] employed a three-layer NN model. Although these nonlinear modeling approaches achieve high predictiv e performance, they pose challenges related to computational cost and training complexity . C . Contribution and no velty This paper presents a bilinear Koopman realization tech- nique from input-output data, incorporating the optimization of the lifting function to enhance its applicability to industrial systems. In the proposed optimization algorithm, RBFs are employed as lifting functions, and their design parameters are optimized using a meta-heuristic approach, specifically PSO. The contrib utions and no velties of this paper are summarized as follo ws: Bilinear K oopman r ealization fr om input-output data: For MIMO nonlinear systems, the prediction performance of an L TI K oopman realization is limited; therefore, a bilinear for- mulation is adopted for Koopman identification. In industrial applications, it is often the case that not all state variables can be measured and only input-output data are av ailable, which is also assumed in this study . By incorporating time-delay co- ordinates into the measured input-output data and defining the input and output delays as embedded states, an input-output bilinear Koopman realization is constructed. For autonomous systems, it is natural that the embedded states are simply lifted since the states include only the output delays. For controlled systems, the selection of the embedded states lifted is important since the states include input delays in addition to output delays. Prior research (e.g., [2], [28]) employs just embedded states with outputs and inputs delay as arguments of lifting functions. Howe ver , the evolution of the inputs is not of interest from the perspectiv e of system identification and feedback control. This is because the inputs are treated as ex- ogenous signals supplied as random or feedback control inputs, rather than as variables e volving autonomously according to a dynamical flow . In this paper , we examine a proper selection of lifting arguments for input-output Koopman realizations. The proposed generalized bilinear K oopman realization solely from input-output data has not been considered in prior works. Optimization of lifting function: The design of lifting func- tions significantly influences the prediction performance of K oopman realization and typically requires substantial manual effort [35]. In this study , RBFs are adopted as the lifting func- tions, and they are optimized using a metaheuristic approach based on multi-step prediction ev aluation. This ev aluation ensures high predictiv e accuracy of the K oopman model. The proposed approach has not been explored in prior re- search and offers se veral advantages: it does not require prior physical insight, leverages PSO as a well-established global optimization solv er , and av oids the inherent challenges of neural network training. Compared to neural network training, PSO offers a robust global optimization framework with easier hyperparameter tuning, clearer objectiv e function definition, and simpler implementation. Application to the real-world industrial system: The pro- posed method is validated through its application to a diesel engine airpath system, which exhibits MIMO characteristics, strong interaction effects, and nonlinear behavior . The method provides a global bilinear formulation that achie ves a balance between predictive accuracy and computational ef ficiency . This study represents the first application of K oopman-based modeling to a diesel engine airpath system, to the best of the 3 authors’ knowledge. Accordingly , the proposed IO-Koopman modeling approach is verified through this real-world indus- trial application. D . Ar ticle organization This paper is org anized as follows. Section II introduces the basic concept of the K oopman operator for input-dependent nonlinear systems, along with methods for representing L TI and bilinear systems using K oopman matrices with finite- dimensional approximations, and techniques for identifying these matrices. Section III proposes a hyperparameter iden- tification method for RBFs, which are employed as lifting functions, in addition to input-output bilinear Koopman mod- eling. Sections IV and V ev aluate the effecti veness of the proposed method through simulation and experimental tests, respectiv ely . The controlled system is the intake and exhaust system of an internal combustion engine, which exhibits strong nonlinearity and MIMO characteristics common in industrial applications. Finally , Section VI summarizes the paper and discusses future research directions. E. Notation The symbol R n denotes the set of real numbers in n - dimensional space, and the symbol R n × m denotes the space of n -row and m -column matrices (or two-dimensional arrays) composed of real numbers. The symbols M and N denote smooth manifolds. The symbol ◦ denotes function composi- tion, the symbol ⊙ denotes the elementwise product, and the symbol ⊗ denotes the Kronecker product. The symbol ∥·∥ F denotes the Frobenius norm defined as ∥ A ∥ F = q P i,j a 2 ij , where a ij is the element of matrix A . The symbol I m denotes the identity matrix of size m × m . The signal sequence is expressed as ( x i ) ∞ i =0 . The space for all signal sequences is expressed as ℓ ( X ) = { ( x i ) ∞ i =0 | x i ∈ X } . The map vec : ℓ ( X ) → R ∞ con vets a signal sequence to vectorized form. The symbol 0 m × n is the zero matrix of size m × n . The symbol † denotes the Moore-Penrose pseudoinv erse of a matrix. The notation ψ [1: n ] is defined as ψ [1: n ] = [ ψ 1 , . . . , ψ n ] ⊤ , where ψ [1: n ] is a vector-v alued function composed of functions ψ i for i ∈ { 1 , 2 , . . . , n } . I I . D A TA - D R I V E N M O D E L I N G V I A K O O P M A N O P E R A T O R A. K oopman operator f or non-autonomous system W e describe the K oopman operator theory for non- autonomous systems [2], [3], [4]. Consider a non-autonomous nonlinear system that includes control inputs: x k +1 = f ( x k , u k ) (1) where f : M × N → M is the mapping, x k ∈ M ⊂ R m , u k ∈ N ⊂ R m denotes the control input at time step k , and N is a smooth manifold. For the system, an extended state is defined as χ k =  x k ν k  (2) where ν k = v ec  ( u i ) ∞ i = k  denotes the state variable con- structed from input sequences ( u i ) ∞ i = k ∈ ℓ ( U ) . The dynamics of the system with the e xtended state are gi ven as χ k +1 = F ( χ k ) =  f ( x k , ν k (0)) S ν k  (3) where ν is the shift operator defined as S ν k = ν k +1 and ν k (0) = u k denotes the first element of nu k at time k . Therefore, the Koopman operator K : F → F for nonlinear systems with inputs is defined as: ( K ψ )( χ k ) = ψ ( F ( χ k )) . (4) The Koopman operator for controlled systems enables a nonlinear system to be represented in a (generally infinite- dimensional) linear form by lifting its state space to an infinite- dimensional space using the observable. The linearity of the K oopman operator helps the use of various linear control theories [2]. Remark 1. The manifolds M and N can be consider ed as M ⊂ R n and N ⊂ R m for most engineering pr oblems [4], [20]. Actually , the states and inputs of the diesel engine airpath system are included in R n and R m , as shown in the simulation and experimental test sections. B. Finite-dimensional K oopman operator In K oopman operator theory , the operator acts on an infinite- dimensional space; therefore, for numerical computations, it is necessary to introduce a finite-dimensional approximation of the Koopman operator . In this section, we present the finite- dimensional approximation of the K oopman operator . The time ev olution of a vector -valued lifting function ψ : R n × R m → R p + q with ψ ∈ F is giv en as ψ ( x k +1 , u k +1 ) = K ψ ( x k , u k ) (5) where K ∈ R ( p + q ) × ( p + q ) is the finite-dimensional Koopman matrix. The lifting function ψ : R n × R m → R p + q , which corresponds to the dictionary of observ ables, can be expressed using two lifting functions ψ x : R n → R p and ψ u : R n × R m → R q as follo ws: ψ ( x k , u k ) =  ψ x ( x k ) ψ u ( x k , u k )  . (6) Since the lifting function ψ x depends on the state variables x , and ψ u depends on both the state variables x and the input u , the generality is not lost [36]. Under this assumption, the time ev olution of the lifting function using the finite-dimensional K oopman operator (i.e., the Koopman matrix K ) is giv en as  ψ x ( x k +1 ) ψ u ( x k +1 , u k +1 )  = K  ψ x ( x k ) ψ u ( x k , u k )  =  K 11 K 12 K 21 K 22   ψ x ( x k ) ψ u ( x k , u k )  (7) where K 11 ∈ R p × p , K 12 ∈ R p × q , K 21 ∈ R q × p , and K 22 ∈ R q × q . Since we are interested in the time e volution of the states, the abov e equation can be simplified as follows: ψ x ( x k +1 ) = K 11 ψ x ( x k ) + K 12 ψ u ( x k , u k ) . (8) 4 Previous studies [18], [19] have demonstrated that various sys- tem representations can be achie ved by appropriately defining the lifting function ψ u ( x k , u k ) . 1) LTI K oopman form: When representing the system in an L TI Koopman form, setting ψ u ( x k , u k ) = u k yields the following L TI-Koopman form: ψ x ( x k +1 ) = K 11 ψ x ( x k ) + K 12 u k . (9) Here, let A = K 11 , B = K 12 , and z k = ψ x ( · ) . Then, the L TI state-space equation for lifted states z k is gi ven as z k +1 = Az k + B u k . (10) 2) Bilinear K oopman form: According to the literature [18], [22], a bilinear K oopman form can be expressed by defin- ing the lifting function ψ u : R n × R m → R m + mp as follo ws: ψ u ( x k , u k ) =  u k ψ x ( x k ) ⊗ u k  (11) where ψ x ( x k ) ⊗ u k is calculated as ψ x ( x k ) ⊗ u k =  u k (1) ψ x ( x k ) , u k (2) ψ x ( x k ) , . . . . . . , u k ( m ) ψ x ( x k )  ⊤ ∈ R mp . (12) Then, the time ev olution of the lifting functions is expressed as   ψ x ( x k +1 ) u k +1 ψ x ( x k +1 ) ⊗ u k +1   =   A B 0 B ∗ ∗ ∗ ∗ ∗ ∗     ψ x ( x k ) u k ψ x ( x k ) ⊗ u k   (13) where A ∈ R p × p , B 0 ∈ R p × m , B =  B 1 , · · · , B m  ∈ R p × mp , and B i ∈ R p × p ( i = 1 , . . . , m ) . Since we are interested in the time ev olution of the states and not in the time ev olution of the inputs, the elements of the K oopman matrix related to the input dynamics can be ignored. Therefore, the abov e equation (13) is reduced to ψ x ( x k +1 ) = A ψ x ( x k ) + B 0 u k + B  ψ x ( x k ) ⊗ u k  . (14) Let z k = ψ x ( · ) denote the lifted state at time k . Then, the bilinear state equation can be written in the follo wing form: z k +1 = Az k + B 0 u k + B ( z k ⊗ u k ) . (15) C . Learning of Koopman matr ix For a nonlinear system with inputs, we assume to measure time-series data of the measurable state variables x k , and the inputs u k are collected, i.e., D = { ( x k , u k ) | k = 1 , . . . , N d } . From the dataset, the snapshot data are set as: Z + and Z w are giv en as Z + =  z 2 , z 3 , . . . , z N d  Z w =   z 1 z 2 · · · z N d − 1 w 1 w 2 · · · w N d − 1 z 1 ⊗ w 1 z 2 ⊗ w 2 · · · z N d − 1 ⊗ w N d − 1   . (16) Then, the state-space matrices A , B 0 , B for the bilinear K oop- man realization are obtained from the following optimization problem: min H   Z + − H Z w   F (17) where H = [ A, B 0 , B ] ∈ R p × p ( m + l +1) . This least-squares problem can be solv ed using the Moore-Penrose pseudoin verse as follo ws: H = Z + [ Z w , U ] † . (18) In the case of a linear system, the elements of the bilinear term are remov ed from the snapshot data. I I I . P RO P O S E D M E T H O D O L O G Y A. Time-dela y coordinates W e explain input-output Koopman modeling [2], which ad- dresses the case where only input-output data can be measured while the system states are unkno wn or unobservable. Since the diesel engine airpath system which has exogenous inputs is considered in this paper, this section handles the nonlinear system with control and exogenous inputs, given as x k +1 = f ( x k , u k , d k ) y k = h ( x k , u k , d k ) (19) where d k ∈ R l is the measurable exogenous input, y k ∈ R n h is the measurable output, and f : R n × R m × R l → R n , h : R n × R m × R l → R n h . For this system with inputs, we assume to obtain the dataset: D = { y k , w k | k = 1 , . . . , N d } , where w k is defined as w k = [ u k d k ] ⊤ . The time-delay coordinates are introduced to address the constraint that only input-output data can be measured. The embedded state at time k , considering delay-step n d , is defined as: ζ k =  ( ζ y k ) ⊤ ( ζ w k ) ⊤  ⊤ ∈ R n ζ (20) with ζ y k =  y ⊤ k y ⊤ k − 1 · · · y ⊤ k − n d  ⊤ ∈ R ( n d +1) n h (21) ζ w k =  w ⊤ k − 1 w ⊤ k − 2 · · · w ⊤ k − n d  ⊤ ∈ R n d ( m + l ) (22) where n ζ = n d ( n h + m + l ) n h . B. IO generalized bilinear K oopman f or m As well as Section II-B, the K oopman model is giv en as ψ x ( ζ k +1 ) = K 11 ψ x ( ζ k ) + K 12 ψ w ( ζ k ) (23) where the lifting function ψ w : R n ζ × R n d ( m + l ) → R ( m + l )+( m + l ) q is defined as: ψ w ( ζ k , w k ) =  w k ϕ w ( ζ k ) ⊗ w k  (24) where ϕ w : R n ζ → R n w is the function. Then, analogous to (13), the standard form is formulated as   ψ x ( ζ k +1 ) w k +1 ϕ w ( ζ k +1 ) ⊗ w k +1   =   A B 0 B ∗ ∗ ∗ ∗ ∗ ∗     ψ x ( ζ k ) w k ϕ w ( ζ k ) ⊗ w k   (25) Since the time ev olution of interest corresponds to the embed- ded states, the abov e equation is simplified as: ψ x ( ζ k +1 ) = Aψ x ( ζ k ) + B 0 w k + B ( ϕ w ( ζ k ) ⊗ w k ) . (26) 5 By defining z w,k = ϕ w ( · ) , the generalized bilinear forms can be re written by z k +1 = Az k + B 0 w k + B ( z w,k ⊗ w k ) . (27) Remark 2. The IO generalized bilinear K oopman form (26) can be vie wed as a bilinear-like LPV Koopman realization by considering ϕ w ( ζ k ) as a scheduling function. C . Argument selection of lifting function for IO-K oopman realization The lifting function in previous studies on the standard IO- K oopman form is defined as ψ x ( ζ k ) =  ζ k ψ [1: n l ] ( ζ k )  ∈ R n ζ + n l (28) where ψ [1: n l ] = [ ψ 1 , . . . , ψ n l ] ⊤ , and ψ i is a scalar function for i ∈ { 1 , . . . , n l } . In previous IO-K oopman realization [2], the lifting function is typically chosen as ψ i ( ζ k ) . Under this lifting function setting, if ψ ( ζ k ) = [ y k , w k − 1 , w k − 1 ⊙ w k − 1 ] ⊤ is giv en, the time ev olution of ψ ( ζ k ) includes inputs. Howe ver , we are not interested in predicting observables that include inputs, since inputs are externally provided—either randomly or through a control law . Thus, employing ψ ( ζ k ) may result in the failure of capturing the essential system dynamic characteristics due to spurious dynamics of exogenous input, which deteriorate the prediction performance or require a significantly larger amount of training data. Although infinite input sequences are assumed in pre vious research, it is not feasible to co ver all possible input sequences due to the finite amount of data in practical situations. Therefore, unlike prior studies [2], [28], we propose explicitly excluding w k from nonlinear lifting to pre vent the identification of spurious dy- namics induced by exogenous inputs. Specifically , we propose the follo wing observable selection: ψ x ( ζ k ) =  ζ k ψ [1: n l ] ( ζ y k )  ∈ R n ζ + n l . (29) Then, the IO bilinear Koopman model is expressed as (30) shown at the bottom of the page. This formulation indicates that the embedded state with input delay is generated by a giv en input. In other words, the input after transitions, w k , is included in the time-e volv ed state z k +1 (= ψ x ( ζ y k +1 )) and the input delay is generated through this input transition. Therefore, the matrices A and B associated with the input- output delay transition are trivially determined by the shift relationship. If there is an observable associated with the input (e.g., y k ⊙ w k − 1 ⊙ w k − 1 ), the observ able is not generated from the giv en input w k . Thus, the matrix A ψ ⋆⋆ in the K oopman matrix A associated with the time ev olution of the lifted state is not uniquely determined. Therefore, in this study , the argument of ψ i is ζ y k instead of ζ k . On the other hand, because ϕ w does not contribute to the time evolution of ψ x ( ζ k ) , we employ a lifting function related to ζ k to increase the flexibility of our representation: ϕ w ( ζ k ) =  ϕ [1: n w ] ( ζ k )  ∈ R n w (31) where ϕ i : R n ζ → R for i ∈ { 1 , · · · , n w } . The output equa- tion is also gi ven as y k = C z k with C =  I n h 0 n h × ( n l − n h )  . W e summarilize the above discussion as follows: Theorem 1. Consider the IO generalized bilinear K oopman r ealization (26) with the state-transition structure in (30), wher e z k (= ψ x ( · )) r epr esents “state” variables updated by a time-in variant linear map A , and the ne w input w k is injected via B 0 and B . Suppose that ψ x contains an input- dependent nonlinear component ψ i ( ζ k ) , such as w k − 1 ⊙ w k − 1 or y k ⊙ w k − 1 . Then the next-time quantity g k +1 (:= ψ i ( ζ k +1 )) corr esponding to ψ i generally depends on the new input w k , and it is structurally impossible to generate g k +1 by the block A ψ ⋆⋆ , and including such nonlinearities in ψ x is causally inconsistent with (30). Pr oof: Assume g k +1 = a ⊤ z k , where a ⊤ is the row of A ψ ⋆⋆ corresponding to the observable ψ i . Fix ζ k (thus z k ) and choose w (1) k  = w (2) k . Since g k +1 depends on w k (e.g., y k +1 ⊙ w k ), one obtains g (1) k +1  = g (2) k +1 , whereas a ⊤ z k is identical for both choices ( ∵ it contains no w k ). This contradiction proves that g k +1 cannot be produced by A ψ ⋆⋆ . Remark 3. Theor em 1 indicates that using ψ i ( ζ k ) may pre vent the essential system dynamics fr om being captur ed due to spurious dynamics induced by the e xogenous input. Corollary 1. T o pr eserve causal consistency with the state- transition structur e in (30), input-dependent nonlinearities must be allocated to ϕ w ( ζ k ) , and the ar guments of ψ i should be r estricted to ζ y k .                 y k +1 y k . . . y k − n d +1 w k w k − 1 . . . w k − n d +1 ψ [1: n l ] ( ζ y k +1 )                 =                  A y y 1 A y y 2 · · · A y yn d A y w 1 A y w 2 · · · A y wn d A y ψ I 0 0 0 0 0 0 0 0 0 . . . 0 . . . 0 0 0 0 0 0 0 I 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 I 0 0 0 0 . . . . . . . . . . . . 0 . . . 0 0 0 0 0 0 0 0 0 I 0 0 A ψ y 1 A ψ y 2 · · · A ψ yn d A ψ w 1 A ψ w 2 · · · A ψ wn d A ψ ψ                                  y k y k − 1 . . . y k − n d w k − 1 w k − 2 . . . w k − n d ψ [1: n l ] ( ζ y k )                 +                 0 0 . . . 0 I 0 . . . 0 0                 w k + B  ϕ w ( ζ k ) ⊗ w k  (30) 6 Pr oof: According to Theorem 1, placing input nonlinear- ities inside ψ x forces their next-time image to depend on w k , which cannot be realized by the state-transient matrix A . Thus, to preserve the interpretation that ψ x represents state variables updated solely by A , input-dependent nonlinearities should be handled through ϕ w ( ζ k ) so that their w k -dependence is routed via B 0 and B , while ψ x is restricted to ζ y k . D . IO generalized K oopman realization This section describes IO generalized bilinear K oopman realization. The snapshot data for the extended states and inputs are set as X =  ζ 1 , ζ 2 , . . . , ζ N d − n d  (32) W =  w 1 , w 2 , . . . , w N d − n d  (33) with ζ i =  ( ζ y i ) ⊤ ( ζ w i ) ⊤  ⊤ ∈ R n ζ (34) ζ y i =  y ⊤ i,n d y ⊤ i,n d − 1 · · · y ⊤ i, 0  ⊤ ∈ R ( n d +1) n l (35) ζ w i =  w ⊤ i,n d − 1 w ⊤ i,n d − 2 · · · w ⊤ i, 0  ⊤ ∈ R n d ( m + l ) . (36) Then, the snapshot data Z + and Z w for IO generalized koopman realization are e xpressed as Z + =  z 2 , z 3 , . . . , z N d  Z w =   z 1 z 2 · · · z N d − 1 w 1 w 2 · · · w N d − 1 z w, 1 ⊗ w 1 z w, 2 ⊗ w 2 · · · z w,N d − 1 ⊗ w N d − 1   (37) and, the IO generalized bilinear Koopman model is obtained from the optimization probelm (17). E. Optimization of lifting functions This section presents the optimization method of lifting functions. Herein, we optimize the design parameters of RBFs. The lifting function is generally selected, such as mechanics- inspired selection, empirical selection (e.g., monomials, poly- nomials), and RBF selection with randomly selected design parameters. Such selections may not provide good perfor- mance since the lifting function is not optimized. The finding method of the proper lifting function is mainly a NN-based approach. This approach has the same difficulties as general NN learning, including hyperparameter setting, ov erfitting, model complexity , computational cost, vanishing and explod- ing gradient problem. In this paper , we adopt a metaheuristic optimization (e.g., PSO) approach, as an alternati ve approach to the NN-based method. The optimization algorithm is shown in Algorithm 1. In the proposed optimization method, RBFs are adopted as the lifting functions. For e xample, polyharmonic RBFs are given as ψ ( x 1 ; θ 1 ) = ∥ x 1 − θ 1 ∥ 2 log ∥ x 1 − θ 1 ∥ 2 (38) ϕ ( x 2 ; θ 2 ) = ∥ x 2 − θ 2 ∥ 2 log ∥ x 2 − θ 2 ∥ 2 (39) where x i and θ i ( i = 1 , 2 ) are the input and center vectors of the RBFs. The center is considered the design parameter . Then, the design parameters of the RBFs are optimized by global optimization such as PSO. Algorithm 1 Optimization of lifting functions 1: Inputs: 2: Max iteration N max 3: RBFs ψ ( x 1 ; θ 1 ) , ϕ ( x 2 ; θ 2 ) 4: Dataset D = { y k , u k , d k | k = 1 , . . . , N d } 5: Outputs: 6: Design parameter vector θ = [ θ 1 , θ 2 ] 7: Koopman matrices A, B 0 , B , C 8: Set initial θ , Maximum iterations N max 9: iteration ← 1 10: while (iteration ≤ N max ) do 11: Compute K oopman matrices A, B 0 , B , C for ψ ( x 1 ; θ 1 ) , ϕ ( x 2 ; θ 2 ) via K oopman realization 12: Update θ by minimizing J ( θ ) 13: if J ( θ ) < ε then 14: break 15: end if 16: iter ation ← iter ation + 1 17: end while 18: Return θ ∗ = arg min J ( θ ) and matrices A, B 0 , B , C As sho wn in Algorithm 1, the proposed optimization process of the centers of RBFs encompasses the realization of the K oopman matrix. The optimization problem is formulated as min θ J ( θ ) , J ( θ ) = N d X k =1 ( y k − ˆ y k ( θ )) 2 (40) where ˆ y k denotes multi-step (not one-step) prediction giv en by ( ˆ z k +1 = A ˆ z k + B 0 w k + B ( ˆ z w,k ⊗ w k ) ˆ y k = C ˆ z k (41) with the initial lifted states z 0 = ψ ( ζ 0 ) and z w, 0 = ϕ ( ζ 0 ) . Matrices A, B 0 , B , C can be obtained by K oopman realization for candidate centers, as sho wn in Section III-D. F . Model v alidation In model validation, the modeling accuracy is ev aluated using the coef ficient of determination, R 2 , of y , defined as R 2 = 1 − P N d k =1 ( y k − ˆ y ( k )) 2 P N d k =1 ( y ( k ) − ¯ y ) 2 (42) where ¯ y is the mean v alue of the N d data [37]. An R 2 score equal to 1.0 indicates that the identified model best fits the target system. From an engineering perspecti ve, an R-squared score is acceptable when it reaches a value from 0.9 to 1.0 [38]. A negati ve value means that the inferred model has a very low ability to represent an equi valent dynamical system. In nonlinear systems, multi-step prediction may not be achiev able even when accurate one-step prediction is achieved [39] . Therefore, we ev aluate both one-step and multi-step predictions. The one-step prediction is defined as ( ˆ x ( t + 1) = f  x ( t ) , u ( t ) , d ( t )  ˆ y ( t ) = h  ˆ x ( t )  (43) 7 Fig. 1. Diesel engine air path system. where ˆ y is the predicted value of y . The one-step ahead is predicted from the giv en input and output data. The multi- step prediction is defined as ( ˆ x ( t + 1) = f  ˆ x ( t ) , u ( t ) , d ( t )  ˆ y ( t ) = h  ˆ x ( t )  (44) with ˆ x (0) = x (0) . The multi-step ahead is predicted from the giv en input and the predicted past output. I V . S I M U L A T I O N A. T arget system The intake and exhaust system of a four-c ylinder diesel engine [35] is shown in Fig. 1. This system includes an ex- haust gas recirculation (EGR) system and a variable geometry turbocharger (VGT). The EGR system regulates the oxygen concentration entering the cylinders by mixing oxygen-lean exhaust gas with fresh air . This reduces the oxygen concentra- tion, lowers the peak comb ustion temperature, and suppresses the formation of harmful nitrogen oxides, which are produced in large quantities at high temperatures. The VGT system adjusts the pressure inside the intake manifold. By narrowing the spacing of the mov able vanes, the flo w path is restricted and the flow rate increases, achieving a supercharging effect ev en at low engine speeds. At high engine speeds, the vane spacing is widened to facilitate exhaust gas flo w . Fresh air is drawn in from the atmosphere and compressed by the compressor . After cooling in the intercooler , it passes through the intake throttle and enters the intake manifold. Exhaust gas also flows into the intake manifold after passing through the EGR valv e. The oxygen concentration in the intake manifold is controlled by operating the EGR valv e. The mixture then flows into the cylinders, and after combustion, the exhaust gas is discharged into the exhaust manifold. The gas leaving the exhaust manifold splits into two paths: one recirculated through the EGR cooler , and the other passing through the turbine and e xiting as e xhaust gas. The burned gas drives the turbine, and the turbine speed is controlled by adjusting the vane angle. In this study , modeling is performed using a data- driv en approach; therefore, further e xplanation of physical modeling is omitted. For details, refer to [40]. B. Simulation setting T raining and test data were generated using a design of experiments (DoE). Refer to [35], [41]. The generated input signals were applied to the plant inputs (i.e., VGT valv e position, EGR valve position, fuel injection amount, and engine speed), and the outputs (i.e., intake manifold pressure T ABLE I S I M U L ATI O N C A S E S O F K O O P M A N R E A L I Z ATI O N . Case Lifting function ψ x Lifting function ϕ w Opt. H-DMD ζ k – – LK  ζ k ψ [1:10] ( ζ k )  – with BLK(poly) 10th-order polynomial of ζ k Similar to ψ x – BLK  ζ k ψ [1:10] ( ζ k )  Similar to ψ x with GBLK  ζ k ψ [1:5] ( ζ y k )   ζ k ϕ [1:5] ( ζ k )  with and EGR rate) were measured. The sampling period was set to 0.1 seconds, resulting in 25,000 data points. The time-series data were normalized as a preprocessing step. T o optimize the RBF hyperparameters, PSO, a reliable global optimization solver , was employed. A delay step of 2 was used. PSO was implemented using the MA TLAB’ s particleswarm function with a maximum of 30 iterations and a swarm size of 300. Parallel computation w as also utilized. T able I shows the koopman forms considered. The number of optimization parameters was set to 10 for all cases. In the table, H-DMD represents Hankel-DMD [42]: z k +1 = Az k + B w w k , LK represents the L TI Koopman realization [2], [14]: z k +1 = Az k + B w w k , BLK represents the standard bilinear K oopman realization [18], [35]: z k +1 = Az k + B w w k + B z ( z k ⊗ w k ) , GBK represents the generalized bilinear Koopman realization (the proposed method): z tk +1 = Az t + B w w t + B z ( z w,k ⊗ w k ) . C . Results and discusstions Fig. 2 illustrates the training and test data generated by applying input signals and measuring the corresponding out- puts. The orange and blue lines represent the training and test datasets, respectively . In the figure, y 1 , y 2 , u 1 , u 2 , d 1 , and d 2 denote the intake manifold pressure [kPa], EGR rate [%], VGT v ane position [%], EGR v alve position [%], fuel injection amount [mm 3 /st], and engine speed [rpm], respectiv ely . Fig. 3 presents the time-series data of prediction performance for test data. The dotted lines represent the predicted v alues, and the solid lines represent actual values. T able II reports the R-squared scores for one-step and long-term predictions for the abov e cases. W e compared our method with con ventional approaches, including Hankel DMD, the L TI K oopman model, the bilinear Koopman model with polynomial basis functions, and the bilinear Koopman model with randomly assigned RBF centers. Among these, Hankel DMD without lifting exhibited the poorest predicti ve performance, and the L TI Koopman model without RBF optimization also performed poorly . Al- though RBF optimization significantly improv ed the predictive accuracy of the L TI Koopman model, the improvement was still insufficient, highlighting the need for identification in the bilinear model. Next, we e xamined the case where the lifting function for the bilinear form w as polynomial and found that its predictive performance was substantially poor . For the standard bilinear form, the R 2 value e xceeded 0.8 after RBF optimization. Howe ver , increasing the number of RBFs led to a marked deterioration in predicti ve performance on the 8 T ABLE II R - S Q UA RE D O F O N E - S T E P A N D L O N G - T E R M P R E D I C TI O NS I N T H E S I M U L ATI O N . Case One-step prediction Long-term prediction learning data learning data test data y 1 y 2 y 1 y 2 y 1 y 2 H-DMD 0.994 0.991 − 14 . 5 − 7 . 85 − 15 . 0 − 7 . 89 LK 0.992 0.989 0.793 0.910 0.778 0.908 BLK(poly) 0.996 0.993 − 14 . 3 − 7 . 77 − 14 . 8 − 7 . 78 BLK 0.996 0.993 0.874 0.988 − 0 . 295 0.347 GBLK 0.995 0.992 0.931 0.951 0.922 0.949 Fig. 2. The learning and test data in simulation. test data, indicating poor generalization. This is likely due to the inclusion of inputs in the observables, as discussed in Section III-D. Finally , the proposed generalized bilinear form B with RBF optimization achiev ed the highest prediction accuracy in this study , with an R 2 value of ov er 0.9, meeting the requirements for engineering applications [38]. Simulation results demonstrated that the proposed method significantly outperformed con ventional approaches. The pre- dicted output of the linearly identified model exhibited poor long-term predicti ve performance in ODE simulations. Diesel engine airpath systems exhibit complex nonlinear character- istics, yet linear models and Hankel DMD rely on linear realizations and therefore fail to deli ver satisfactory results. W e then examined the bilinear K oopman realization. First, for the bilinear K oopman model using polynomial basis functions– similar to previous studies [24], [25]–the lifting functions were selected based on findings from SINDy modeling [35]. Howe ver , this approach prov ed ineffecti ve for modeling diesel engine airpath systems. Since designers cannot successfully reuse the basis functions employed in SINDy , alternativ e approaches are required. Pre vious literature has used randomly centered RBFs as lifting functions, but as noted in [43], K oop- man models with non-optimized RBFs tend to div erge. Our simulations confirmed this observ ation: while one-step pre- dictions achiev ed very high R-squared exceeding 0.98, long- term predictions failed to meet performance requirements. Our proposed method addresses this issue by optimizing the centers of the RBFs used as lifting functions based on multi-step prediction performance. The bilinear Koopman model with RBF optimization demonstrated superior predictiv e accuracy compared to other methods. Furthermore, we confirmed that prediction and generalization performance depend strongly on the choice of arguments for the lifting function. Pre vious studies typically used embedded state variables consisting of Fig. 3. The long-ter m predictive performance of y 1 and y 2 for test data in simulation. time delays of inputs and outputs, but we found that these methods suffered from poor prediction and generalization. In contrast, our proposed argument selection strategy achiev ed significantly better performance. These results demonstrate the effecti veness of the proposed method. V . R E A L - W O R L D E X P E R I M E N T A L V E R I FI C A T I O N A. Experimental setting A real-world engine bench test, as shown in Fig. 4, was conducted. The experimental environment is the same as in the literature [44]. The v ariable y 2 represents the mass air flow (MAF) [mg/st], which is related to the EGR ratio, while the other variables are the same as those used in the simulation. In the engine bench tests, training and test data are collected under closed-loop conditions. The sampling period was 100 ms. In contrast, the simulation section used open- loop data, which are generally considered preferable. Here, we demonstrate that the proposed method remains effecti ve in closed-loop experiments, where measurement noise tends to be amplified. The controller employed in the closed-loop experiments w as a PID-based mass-production controller . B. Results and discussions Fig. 5 shows the the learning and test data. The K oopman models were deri ved from the learning data. Fig. 6 sho ws the Fig. 4. Engine bench test. Fuel type:diesel, Arrangement cylinders: 4 in line, Maximum power:110 kW (150 PS)3600 r pm, Maximum torque: 350 Nm (35.69 kgm)1800-2600 r pm (net), Combustion type: Direct injection with water-cooled 4-valv e DOHC (double overhead camshaft). 9 T ABLE III R - S Q UA RE D O F O N E - S T E P A N D L O N G - T E R M P R E D I C TI O NS I N T H E E X P E R I M EN TAL T E S T . Case One-step prediction Long-term prediction learning data learning data test data y 1 y 2 y 1 y 2 y 1 y 2 LK 0.999 0.994 0.988 0.973 0.917 0.834 BKL 0.999 0.997 0.997 0.992 − 0 . 774 − 3 . 68 GBKL 0.999 0.997 0.996 0.992 0.936 0.911 Fig. 5. The learning and test data in experimental test. prediction performance of the proposed method on the test data. These signals are normalized. In addition, we compare the performance of the conv entional and proposed methods. T able III reports the R-squared values for one-step predictions on the learning data and for both one-step and long-term predictions on the test data. The table also compares the results of the linear and bilinear K oopman realizations. The cases presented in T able III corresponds to those shown in T able I. From the results, the proposed method provides the best performance among them, as well as simulation verification. The L TI K oopman form is limited to the pre- diction performance. The predictions of the standard bilinear K oopman form, where lifting functions are ψ x = ψ w div erge for the test data. Ho wev er, the proposed generalized bilinear K oopman form pro vides the desired performance for learning and test data. From the above, the ef fectiveness of the proposed methodology was e xperimentally verified. Section V has examined the effecti veness of the proposed method through an actual engine bench test. Similar to the simulation verification, the proposed method–namely , the IO bilinear K oopman realization–achiev es the best performance among the compared methods, including the L TI K oopman re- alization and the standard bilinear Koopman realization. In the simulation section, the learning data were acquired through an open-loop test. General system identification typically prefers rich excitation signals under open-loop testing conditions. In the experimental section, we assumed a more realistic scenario and obtained transient data from closed-loop experiments as training data. In other words, we considered that in actual tests, applying random vibrations using open-loop tests would not be easy and would require additional technical expertise. Even under such conditions, we demonstrated that the proposed method can achieve the desired results for both the training and test data. Fig. 6. The long-ter m predictive performance of y 1 and y 2 for test data in experimental test. V I . C O N C L U S I O N S This study presented an input-output bilinear Koopman realization with a lifting function optimization algorithm. The proposed framework integrates bilinear Koopman realization, input-output data utilization, and lifting function optimization. This approach enables K oopman-based modeling under prac- tical constraints where only input-output data are av ailable, significantly expanding its applicability to industrial systems. The key contributions of this work include the in vestigation of bilinear form structures for K oopman realization, the system- atic selection of lifting function ar guments when using input- output data, the optimization of lifting functions via PSO, and the successful application of the proposed method to an industrial diesel engine airpath system. Unlike pre vious studies that simply combined embedded states with input and output delays–thus including the time ev olution of uninformative inputs–our method addresses argument selection to enhance predictiv e performance. Furthermore, because the design of the lifting function strongly influences accuracy and is difficult for users to determine, the proposed PSO-based optimization– where the centers of the RBFs employed as the lifting function are optimized–provides a practical and effecti ve solution. The method was validated on a diesel engine intake and exhaust system exhibiting strong nonlinearity and MIMO character- istics. Simulation and experimental test results demonstrated substantial improvements in predictiv e accuracy compared to con ventional approaches, achieving an R-squared over 0.90, which meets engineering requirements. Future work will focus on constructing controller with IO K oopman model. R E F E R E N C E S [1] P . V an Overschee and B. De Moor, “N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems, ” Automatica , v ol. 30, no. 1, pp. 75–93, Jan. 1994. [2] M. K orda and I. Mezi ´ c, “Linear predictors for nonlinear dynamical systems: K oopman operator meets model predictiv e control, ” Automatica , v ol. 93, pp. 149–160, July 2018. 10 [3] M. O. Williams, C. W . Rowley , and I. G. Ke vrekidis, “A kernel- based method for data-driven koopman spectral analysis, ” Journal of Computational Dynamics , vol. 2, no. 2, pp. 247–265, 2015. [4] J. L. Proctor, S. L. Brunton, and J. N. Kutz, “Generalizing K oopman Theory to Allow for Inputs and Control, ” SIAM J ournal on Applied Dynamical Systems , vol. 17, no. 1, pp. 909–930, Jan. 2018. [5] K. Huang, K. W ei, F . Li, C. Y ang, and W . Gui, “LSTM-MPC: A Deep Learning Based Predictive Control Method for Multimode Process Control, ” IEEE Tr ansactions on Industrial Electr onics , vol. 70, no. 11, pp. 11 544–11 554, Nov . 2023. [6] K. Zhang, X. Luan, F . Ding, and F . Liu, “A Bayesian transfer sparse identification method for nonlinear ARX systems, ” International Journal of Adaptive Control and Signal Processing , vol. 38, no. 10, pp. 3484–3502, Oct. 2024, publisher: W iley . [7] S. L. Brunton, J. L. Proctor , and J. N. Kutz, “Sparse Identification of Nonlinear Dynamics with Control (SIND Yc), ” May 2016, arXiv:1605.06682 [math]. [8] D. Bhattacharya, R. Hashem, L. K. Cheng, and W . Xu, “Nonlinear Model Predictiv e Control of a Robotic Soft Esophagus, ” IEEE T ransactions on Industrial Electronics , vol. 69, no. 10, pp. 10 363– 10 373, Oct. 2022. [9] S. L. Brunton, M. Budi ˇ si ´ c, E. Kaiser, and J. N. Kutz, “Modern K oopman Theory for Dynamical Systems, ” SIAM Review , vol. 64, no. 2, pp. 229–340, May 2022. [10] P . Bev anda, S. Sosno wski, and S. Hirche, “Koopman Operator Dynamical Models: Learning, Analysis and Control, ” Annual Reviews in Contr ol , vol. 52, pp. 197–212, 2021. [11] P . J. Schmid, “Dynamic mode decomposition of numerical and experimental data, ” Journal of Fluid Mechanics , vol. 656, pp. 5–28, Aug. 2010. [12] J. L. Proctor , S. L. Brunton, and J. N. Kutz, “Dynamic Mode Decomposition with Control, ” SIAM Journal on Applied Dynamical Systems , vol. 15, no. 1, pp. 142–161, Jan. 2016. [13] M. O. W illiams, I. G. Ke vrekidis, and C. W . Rowley , “A Data–Driv en Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition, ” Journal of Nonlinear Science , vol. 25, no. 6, pp. 1307–1346, Dec. 2015. [14] M. O. Williams, M. S. Hemati, S. T . Dawson, I. G. Ke vrekidis, and C. W . Ro wley , “Extending Data-Driven Koopman Analysis to Actuated Systems, ” IF AC-P apersOnLine , vol. 49, no. 18, pp. 704–709, 2016. [15] Y . Meng, L. Li, X. W ang, W .-W . Huang, E. Q. W u, and L. Zhu, “Online Koopman Operator-Based Feedforward Compensation Strategy for Fast T ool Servos W ith Robust High-Bandwidth Control, ” IEEE T ransactions on Industrial Electr onics , v ol. 72, no. 3, pp. 2958–2967, Mar . 2025. [16] D. Zhao, X. Y ang, Y . Li, L. Xu, J. She, and S. Y an, “A Kalman-K oopman LQR Control Approach to Robotic Systems, ” IEEE T ransactions on Industrial Electr onics , vol. 71, no. 12, pp. 16 047–16 056, Dec. 2024. [17] S. Y u, C. Shen, and T . Ersal, “Autonomous Driving using Linear Model Predictiv e Control with a Koopman Operator based Bilinear V ehicle Model, ” IF AC-P apersOnLine , vol. 55, no. 24, pp. 254–259, 2022. [18] D. Bruder, X. Fu, and R. V asude van, “Advantages of Bilinear K oopman Realizations for the Modeling and Control of Systems W ith Unknown Dynamics, ” IEEE Robotics and A utomation Letter s , v ol. 6, no. 3, pp. 4369–4376, July 2021. [19] D. Goswami and D. A. Paley , “Global bilinearization and controllability of control-affine nonlinear systems: A Koopman spectral approach, ” in 2017 IEEE 56th Annual Confer ence on Decision and Contr ol (CDC) . Melbourne, Australia: IEEE, Dec. 2017, pp. 6107–6112. [20] M. W ang, X. Lou, and B. Cui, “Deep bilinear Koopman realization for dynamics modeling and predictiv e control, ” International Journal of Machine Learning and Cybernetics , v ol. 15, no. 8, pp. 3327–3352, Aug. 2024. [21] D. Zhao, B. Li, F . Lu, J. She, and S. Y an, “Deep Bilinear K oopman Model Predictive Control for Nonlinear Dynamical Systems, ” IEEE T ransactions on Industrial Electr onics , vol. 71, no. 12, pp. 16 077–16 086, Dec. 2024. [22] M. Abtahi, F . M. Araghi, N. Mojahed, and S. Nazari, “Deep Bilinear K oopman Model for Real-Time V ehicle Control in Frenet Frame, ” July 2025, arXiv:2507.12578 [eess]. [23] L. Shi, Z. Liu, and K. Karydis, “Koopman Operators for Modeling and Control of Soft Robotics, ” Current Robotics Reports , vol. 4, no. 2, pp. 23–31, July 2023. [24] X. Zhang, M. Han, and X. Y in, “Reduced-order Koopman modeling and predictive control of nonlinear processes, ” Computers & Chemical Engineering , vol. 179, p. 108440, Nov . 2023. [25] J. W ang, B. Xu, J. Lai, Y . W ang, C. Hu, H. Li, and A. Song, “ An Improved Koopman-MPC Framework for Data-Driven Modeling and Control of Soft Actuators, ” IEEE Robotics and Automation Letters , vol. 8, no. 2, pp. 616–623, Feb . 2023. [26] B. Lusch, J. N. Kutz, and S. L. Brunton, “Deep learning for univ ersal linear embeddings of nonlinear dynamics, ” Nature Communications , vol. 9, no. 1, p. 4950, Nov . 2018. [27] Y . Y ang, Y . Jiang, W . Cao, C. Li, W . Jiang, H. Mo, X. Li, X. Xu, and X. Zhang, “Adaptive Deep Koopman Operators for Soft Robot Control: Application to Luban Lock Disassembly, ” IEEE T ransactions on Industrial Electr onics , pp. 1–12, 2026. [28] L. Xiawen, M. Chetan, C. Shichao, W . Y ajun, and D. L. R. Jaime, “Parameters determination of time-delayed embedding with application to K oopman operator-based model predictiv e frequency control, ” CSEE Journal of P ower and Energy Systems , 2021. [29] R. Moriyasu, S. Nojiri, A. Matsunaga, T . Nakamura, and T . Jimbo, “Diesel engine air path control based on neural approximation of nonlinear MPC, ” Control Engineering Practice , vol. 91, p. 104114, Oct. 2019. [30] D. Zhao, C. Liu, R. Stobart, J. Deng, E. W inward, and G. Dong, “An Explicit Model Predictive Control Framework for T urbocharged Diesel Engines, ” IEEE T ransactions on Industrial Electronics , vol. 61, no. 7, pp. 3540–3552, July 2014. [31] H. Huang, X. Peng, W . Du, S. X. Ding, and W . Zhong, “Nitrogen Oxides Concentration Estimation of Diesel Engines Based on a Sparse Nonstationary Trigonometric Gaussian Process Regression W ith Maximizing the Composite Likelihood, ” IEEE T ransactions on Industrial Electr onics , vol. 70, no. 11, pp. 11 744–11 753, Nov . 2023. [32] X. W ei and L. Del Re, “Gain Scheduled $H {\ infty } $ Control for Air Path Systems of Diesel Engines Using LPV T echniques, ” IEEE T ransactions on Contr ol Systems T echnology , vol. 15, no. 3, pp. 406–415, May 2007. [33] M. Iwadare, M. Ueno, and S. Adachi, “Multi-V ariable Air-Path Management for a Clean Diesel Engine Using Model Predicti ve Control, ” SAE International Journal of Engines , vol. 02, no. 1, pp. 764–773, Apr . 2009. [34] D. Gagliardi, T . Othsuka, and L. D. Re, “Direct C/GMRES Control of The Air Path of a Diesel Engine, ” IF A C Pr oceedings V olumes , vol. 47, no. 3, pp. 3000–3005, 2014. [35] S. Y ahagi, H. Seto, A. Y onezawa, and I. Kajiwara, “Sparse Identification and Nonlinear Model Predictive Control for Diesel Engine Air Path System, ” International Journal of Contr ol, Automation and Systems , vol. 23, no. 2, pp. 620–629, Feb . 2025. [36] S. H. Son, A. Narasingam, and J. S.-I. Kwon, “Handling plant-model mismatch in Koopman L yapunov-based model predictiv e control via offset-free control framework, ” Oct. 2020, arXi v:2010.07239 [cs, eess]. [37] R. Haber and H. Unbehauen, “Structure identification of nonlinear dynamic systems—A surv ey on input/output approaches, ” Automatica , vol. 26, no. 4, pp. 651–677, July 1990. [38] B. Schaible, Hong Xie, and Y ung-Cheng Lee, “Fuzzy logic models for ranking process effects, ” IEEE T ransactions on Fuzzy Systems , vol. 5, no. 4, pp. 545–556, Nov . 1997. [39] Y . Xiao, X. Zhang, X. Xu, X. Liu, and J. Liu, “Deep Neural Networks W ith Koopman Operators for Modeling and Control of Autonomous V ehicles, ” IEEE Tr ansactions on Intelligent V ehicles , vol. 8, no. 1, pp. 135–146, Jan. 2023. [40] S. Kaneko, Y . Y amasaki, H. Ohmori, H. Mitsuo, I. Mizumoto, M. Ichiyanagi, A. Matsunaga, and T . Jimbo, Model Based Contr ol for Automotive engines . CORON A, 2019. [41] A. Y onezawa, H. Y onezawa, S. Y ahagi, I. Kajiwara, S. Kijimoto, H. T aniuchi, and K. Murakami, “Sparse Identification of Nonlinear Dynamics W ith Library Optimization Mechanism: Recursive Long-T erm Prediction Perspectiv e, ” IEEE T ransactions on Cybernetics , 2026. [42] C. C. Zawacki and E. H. Abed, “Dynamic Mode Decomposition for Control Systems with Input Delays, ” IF AC-P apersOnLine , vol. 56, no. 3, pp. 97–102, 2023. [43] H. Shi and M. Q.-H. Meng, “Deep Koopman Operator W ith Control for Nonlinear Systems, ” IEEE Robotics and Automation Letters , vol. 7, no. 3, pp. 7700–7707, July 2022. [44] S. Ishizuka, I. Kajiwara, J. Sato, Y . Hanamura, and S. Hanaw a, “Model-free adaptiv e control scheme for EGR/VNT control of a diesel engine using the simultaneous perturbation stochastic approximation, ” T ransactions of the Institute of Measurement and Contr ol , vol. 39, no. 1, pp. 114–128, Jan. 2017.