Linear model predictive control with lifted bilinear models by Koopman-based approach

This study proposes a linear Model Predictive Control (MPC) method that combines high prediction accuracy with low computational cost, using a lifted bilinear model based on Koopman theory. In recent years, there has been a growing interest in approaches to learning prediction models that determine the performance of MPC through lifting linearization and lifting bilinearization based on Koopman theory. In these methods, a linear or bilinear model reflecting nonlinear characteristics of the target system can be obtained by lifting the states to a higher dimensional space with observable functions (observables). In particular, lifting bilinearization provides more accurate models for the nonlinear input affine system, but when combined with MPC, it requires solving a nonlinear optimization problem, which is computationally expensive. Therefore, in this study, we formulate a linear MPC using a lifted bilinear model that can make accurate predictions for the input affine system, thereby realizing an MPC algorithm with high accuracy and low computational cost. In the proposed method, we first formulate a prediction error correction method for the lifted bilinear model by introducing constraints based on the observables and corrective inputs. Furthermore, we propose a linear MPC that can make accurate predictions using the lifted bilinear model by utilizing prior-predicted states based on the optimal solution at the previous time. We evaluate the effectiveness of the proposed method through numerical simulations using a planar drone.


Introduction
Model Predictive Control (MPC) is gaining attention as a method that enables the systematic design of optimal control controllers for a variety of systems, taking into account constraints on states and inputs [1,2]. Since MPC uses a prediction model of the target system to compute a sequence of control inputs that optimizes the states for a finite time horizon, the prediction accuracy of the states is crucial to achieve high-performance control that can accommodate the objective function.
In recent years, techniques using lifting linearization based on Koopman theory have been actively studied to obtain linear prediction models that reflect the characteristics of nonlinear systems [3]. Lifting linearization is a method that allows the state to be lifted to a higher dimensional space using the nonlinear outputs of the system, called observables, and to represent the time evolution of the observables as linear dynamics. The lifted linear model obtained in this way allows various theories for linear systems that have already been established to be applied to nonlinear systems.
In lifting linearization, a linear model is learned in a data-driven manner based on observables. At this time, it is important to select appropriate observables, and a variety of methods have been proposed for this purpose. One popular method is Extended Dynamic Mode Decomposition (EDMD) [4], which uses a dictionary of generic observables such as monomials up to a certain degree. Also, techniques that optimize or learn the observables itself have been proposed [5,6]. In another direction, there are approaches such as Dual Faceted Linearization (DFL) [7][8][9] that utilize the physical knowledge of the target system to determine the observables.
Here, the challenge was that model learning by lifting linearization may not be sufficiently accurate, especially for input affine systems that describe dynamics such as aerial and ground vehicles. The main reason for this is that these systems add inputs via statedependent rotational matrices, and linear models with state-independent input coefficient matrices cannot accurately represent their characteristics.
To address this problem, approaches have been proposed in which lifted bilinear models are learned for the input affine system in a data-driven manner, leading to nonlinear MPC with performance comparable to that of the case with a full knowledge of the target system [10]. This method exploits the fact that the Koopman Canonical Transformation (KCT) allows the input affine system to be globally transformed into a bilinear system in a high-dimensional space under certain conditions [11]. However, the use of the lifted bilinear models as prediction models for MPC has faced challenges of computational cost because it requires to solve nonlinear optimization problems.
Therefore, in this study, we establish an accurate and low-computational-cost linear MPC algorithm utilizing lifted bilinear models that can make accurate predictions for an input affine system. In the proposed method, we first construct a prediction error correcting MPC that compensates for the prediction error of the lifted bilinear model by introducing constraints on the observables and corrective inputs. Then, we approximate the lifted bilinear model as a time-varying linear model using the prior-predicted states and linearize the nonlinear constraints. The prior-predicted states can be easily obtained by taking the advantage of the MPC algorithm's feature of iteratively solving the optimization on the prediction horizon at each time.
The structure of the remainder of this paper is as follows. Section 2 describes preliminary knowledge of the lifting linearization and bilinearization, and how to identify the lifted bilinear model. In Section 3, we propose a prediction error correction method for MPC using the lifted bilinear model. In Section 4, we extend the method of Section 3 and propose an algorithm for linear MPC utilizing the lifted bilinear model. Section 5 presents the results of a numerical simulation that deals with a planar drone as an example of an input affine system to evaluate the effectiveness of the proposed method. Finally, Section 6 provides conclusion.

Koopman operator
The Koopman operator [12] was originally proposed for autonomous systems as a method to map nonlinear dynamics to an infinite-dimensional linear space. Later, efficient algorithms such as EDMD [4] were proposed to approximate nonlinear dynamics to finitedimensional linear spaces in a data-driven manner.
Consider the discrete-time nonlinear autonomous system where k is the time index and x k ∈ R n . Let ψ : R n → C be observables, which are nonlinear functions of the state. The Koopman operator K is an infinitedimensional linear operator acting on the observables, and the time evolution of the observables according to (1) can be expressed as linear dynamics: Also, for a continuous-time nonlinear autonomous systemẋ with x ∈ R n as state, there exists a Koopman operator semigroup K t corresponding to a sample time t > 0.
The infinitesimal generatorK of this semigroup is which is referred to as the continuous-time Koopman operator. In addition, Koopman eigenfunctions φ : R n → C are functions corresponding to Koopman eigenvalues λ ∈ C ofK that evolves linearly along the solution of (3), whose future time values are represented by the following relation: Although Koopman operator is generally infinite dimensional, it can be approximated as a finite dimensional matrix by data-driven methods such as EDMD using a finite number of selected observables. Such methods have also been extended to be applied to control systems with inputs [13,14]. This allows us to identify the state evolution of a nonlinear system lifted in higher dimensional space by observables as linear dynamics, which is called lifting linearization. Furthermore, various linear control theories can be applied to the obtained linear model.

Lifting bilinearization
In this paper, we consider the following input affine systems as target systems: where x ∈ R n and u ∈ R m . Input affine systems can describe the dynamics of aerial and ground vehicles and mechanical systems. In such systems, the input generally act on state via a state-dependent matrix g(x). For example, in vehicle dynamics, the driving force affects the state via an attitude angle-dependent rotation matrix. When the target system (6) is identified as a linear model using the lifting linearization, the coefficient matrix of the input is obtained as a stateindependent constant matrix. Therefore, the effect of state-dependent coefficients on the input may not be appropriately captured. Especially for systems with large attitude angle changes, it may not be possible to obtain a sufficiently accurate model for the original system.
On the other hand, it is known that input affine systems can be transformed into bilinear systems under certain conditions using the Koopman Canonical Transformation (KCT) [10,11].
Let (λ i , φ i ), i = 1, . . . , n z be eigenvalue-eigenfunc tion pairs of Koopman operator for the autonomous dynamicsẋ = f (x) of (6). Also, let φ i : R n → R and z = T(x) = [φ 1 (x) · · · φ N (x)] T , then KCT is defined as where F ∈ R n z ×n z is a diagonal matrix with λ i in (i, i) elements. Note that in (7a), the state x of the original system (6) is represented by a linear combination of Koopman eigenfunctions. At this time, if the space spanned by the Koopman eigenfunctions φ i , i = 1, . . . , n z forms an invariant subspace of L g i , i = 1, . . . , m, then the input affine system (6) is bilinearizable in n z -dimensional space. As a result, L g i T(x) = G i z, G i ∈ R n z ×n z , and the input affine system can be expressed in the following Koopman Bilinear Form (KBF) (see [11] for details):

Identification of lifted bilinear model
Now we describe how to identify the lifted bilinear model represented by KBF when the target system (6) is unknown. First, assume that the following state and input sequence are obtained as the observation data of the target system: where the sampling interval is t. In the data acquiring experiment, a controller designed for stabilization or tracking may be applied to obtain the data. The data D may be a collection of multiple data sets acquired through separate experiments. In such a case, when constructing matrices Z + , Z u discussed below, care should be taken not to associate data across different data sets. Next, let the vector of observables be h(x) ∈ R n z . Here, the observables may be selected based on any knowledge of the physical model of the target system. In a more general way, a generic basis function, such as a monomial up to a certain degree, may be chosen. In that case, however, a sufficiently high order may have to be chosen to obtain a highly accurate model.
Here, a simple way to satisfy (7a) is to include the original state in the observables. Let η(x) ∈ R n η be a vector of observables other than the original state x, and assume the vector of observables in the form Now we consider the following discrete-time lifted bilinear model with z k as the lifted states: where u k (i) is the ith component of u k . Given data D and calculating z k = h(x k ), k = 0, . . . , M, the coefficient matrices A, B 1 , · · · , B m of (11a) can be obtained by solving the optimization problem (12) (12) is the optimization problem of determining the coefficient matrix so as to minimize the prediction error between samplings. The solution can be obtained by the least-squares method. Let H = [A B 1 · · · B m ] ∈ R n z ×n z (m+1) and construct matrices as Note that each column of Z + and Z u is associated with the left and right sides of (11a) at each sampling. Then, the estimation of H can be calculated bŷ Also, from (10), the matrix C in (11b) can be obtained as

MPC with lifted bilinear model
In this section, we consider MPC using the lifted bilinear model obtained as the prediction model. First, we consider the formulation as a nonlinear MPC and construct the optimization problem to be solved at each time as follows: where N p is the prediction horizon and x t is the observed states at time t. A c , b c are a matrix and a vector constituting linear constraints on the state and the input. The objective function J 1 is given by where S, Q ∈ R n×n , R ∈ R m×m area positive semidefinite weight matrices and x r k , k = 0, . . . N p and u r k , k = 0, . . . , N p − 1 are the target state and the target input, respectively. We use the notation of weighted norm x Q = x T Qx.
In the optimization problem (16), the observed state x t at the current time is lifted to the lifted state space by the vector of observables h(·), and the evolution of the lifted states on the prediction horizon is predicted by the lifted bilinear model. Furthermore, by projecting the predicted lifted states at each step onto the original state space and evaluating it in the objective function, the optimal input sequence u k , k = 0, . . . ., N p − 1 is obtained.

Prediction error correction
The lifted bilinear model is expected to have higher prediction accuracy for the input affine system than the lifted linear model. However, since it is still an approximated model based on the data, some amount of error is inevitable. Furthermore, in this study, the lifted bilinear model is further approximated as a time-varying linear model to formulate a linear MPC, as described in detail in the next section. This may further increase the prediction error.
When the lifted model is used as a prediction model, especially for the nonlinear component η(x) of the original states, the prediction error is likely to be large if the choice of the observables is not adequate. And since the prediction of the original states is calculated using the predictions of the nonlinear elements, the prediction errors may affect them in various ways. As a result, the objective function may not be properly evaluated and the performance of the closed-loop system with MPC may deteriorate.
To compensate for the effects of such errors, we consider state constraints based on the relationship of the observables to ensure that the values of the lifted states are properly balanced on the prediction horizon. Also, to avoid inconsistency in the state-input relationship in the state equation, a corrective input term is added to the component of the lifted states that correspond to nonlinear functions. The optimization problem with these considerations is formulated as follows: where u k ∈ R n η is the corrective input vector for the n + 1 to n z th states at each prediction step k. z k (i : j) is the vector of the i to jth elements of z k following the MATLAB notation. Defining a weight matrix regarding u k as R ∈ R n η ×n η , the objective function J 2 is given by In the optimization problem (17), by (17e), we consider that z k (1 : n) corresponding to the original state and z k (n + 1 : n z ) corresponding to the portion of the lifted states augumented by the nonlinear function vector η(x) satisfy the relationship of the observables at each time step. By solving this optimization problem, the optimal corrective inputs are computed so that the elements of the lifted states at each step satisfy the constraints, by which the prediction error due to the lifted model inaccuracy can be compensated. Since an improvement in the prediction accuracy of z k leads to an improvement in that of x k , the control performance of the MPC is expected to be improved.
Here, the corrective input u k is computed to be as small as possible thanks to the objective function, so the optimization of the state trajectory is realized using the original input as much as possible. In addition, if the lifted model accurately approximates the nonlinear model, the corrective input become small enough. In other words, it is expected that the corrective input u k compensates for the uncertainty in the dynamics of the prediction model as necessary, but does not interfere with the optimization calculation of the original input.
Note that (17c) is a prediction model based on the lifted bilinear model, which encompasses the lifted linear model. Indeed, if the components of η(x) contain a constant and all the corresponding elements of B i other than the constant are 0, then B i z k =B i and the bilinear terms disappear. In such a case, the model becomes linear with respect to the input. Therefore, the proposed prediction error correction MPC can be equally applied to the case where the prediction model is a lifted linear model.

Linearization with prior-predicted states
The problem (17) is a nonlinear optimization problem because it involves nonlinear constraints (17e) in addition to the bilinear prediction model (17c).
To formulate a linear MPC, we consider linearizing a bilinear prediction model (17c) and nonlinear constraints (17e) with the a previously obtained state sequence prediction for the prediction horizon. We call this prior-predicted states. Fortunately, MPC calculates the optimal state sequence at the same time it finds the optimal input sequence at each time. Therefore, the prior-predicted states can be easily obtained by utilizing the optimal state sequence calculated at the previous time. The details of this procedure are described in the next subsection.
Let z * k , k = 0, . . . , N p be the prior-predicted states at the kth step of the prediction horizon. Also, the priorpredicted value of the original states corresponding to z * k can be obtained by x * k = Cz * k , k = 0, . . . , N p . First, the lifted bilinear model for the prediction can be regarded as a time-varying linear system on the prediction horizon, using the prior-predicted states as follows: whereB(z * k ) = [B 1 z * k · · · . B m z * k ]. Next, the nonlinear constraint (17e) on the lifted state at the kth step of the prediction horizon can be reformulated as the following linear constraint by linear approximation around the prior-predicted states: (19) From the above, given the prior-predicted states z * k , k = 0, . . . N p , the optimization problem for the proposed LMPC is formulated as In the optimization problem (20), the only nonlinear calculation is calculation of the lifted states z 0 at the initial step (20b) and the other calculations all consist of a quadratic form objective function and linear constraints. Therefore, the formulation can be converted to the dense form for the quadratic programming problem, and solved with an efficient QP solver.

Obtaining prior-predicted states
Next, we explain how to obtain the prior-predicted states z * k . Given the optimal input sequence u k and u k obtained by solving the problem (20) at a certain time t, using the bilinear prediction model (20d), it is possible to calculate the optimal value of the lifted states at each step z * − k , k = 0, . . . , N p . Provided that the control interval coincides with the sampling interval t of the discrete-time prediction model (20d), the priorpredicted states at the next time t + t can be approximately predicted by shifting the optimal state sequence z * − k at the previous time by one step, except for the N p th step. For the N p th step, the value of (N p − 1)th step can be used. That is, the prior predicted state z * k can be obtained as follows: At the initial time of the control loop, the above method cannot be applied because no previous solution is available. Instead, we can use the lifted state of the system in all steps of the prediction horizon. That is, given the initial state x t0 , the prior-predicted states are given by The prior-predicted states at the initial time also may be obtained by solving a nonlinear optimization problem with a nonlinear model in advance. This is expected to improve the control performance around the initial time because the approximation of the lifted bilinear model and the nonlinear constraints are more accurate. For cases where the control interval and the sample time in the discrete-time prediction model do not match, the optimal state sequence at the previous time can be interpolated to virtually shift the optimal state sequence by the control interval.

Algorithm
From the discussion so far, the proposed linear MPC algorithm based on the lifted bilinear model can be summarized as Algorithm 1.

Algorithm 1 LMPC based on lifted bilinear model
Require: Reference trajectory (x r , u r ) while Controller is running do if t=0 then Set the prior-predicted states z * k by (22) end if Observe current state x t of the system Solve the optimization problem (20) with the prior-predicted states z * k Update the prior-predicted states z * k by (21) Apply first input u 0 to the system end while The major advantage of this algorithm is that it can simultaneously compensate for the prediction error and improve the model accuracy of the approximated lifted bilinear model. In other words, even if there is an error in the prior-predicted states at a certain time and the time-varying linear prediction model is not accurate enough, the resulting prediction error can be reduced by the corrective input so that the constraints on the lifted state are satisfied. Then, by utilizing the optimal state sequence with reduced prediction errors, the approximation accuracy of the lifted bilinear model and constraints on the lifted state is improved. This iterative accuracy improvement effect allows the proposed LMPC algorithm to achieve prediction performance equivalent to that of a lifted bilinear model (11), while being in the form of a quadratic programming problem which can be solved efficiently.

Target system
In this section, we evaluate the effectiveness of the proposed method through numerical simulations. As a target system, we consider a drone dynamics operating in the two-dimensional plane shown in Figure 1.
The system is an input affine nonlinear system whose dynamics is given by  d dt  Table 1.
In the following sections, we first identify a lifted bilinear model for the target system (23). Then, using the obtained model, we apply the LMPC algorithm of the proposed method and evaluate its performance. Two control methods are considered for comparison: nonlinear MPC which uses the true system as the prediction model (NMPC), and linear MPC which uses a linear approximation model around the equilibrium point as the prediction model (Nominal LMPC).
Here, the linear approximation model for the Nominal LMPC is obtained from the first-order Taylor expansion around the equilibrium statex = [0 0 0 0 0 0] T and the equilibrium inputū = [ 1 2 mg 1 2 mg] T , as where

Model identification
To identify the lifted bilinear model, since the target system (23) is unstable, a stabilizing controller to the origin was applied to obtain the observation data. As a stabilizing controller, LQR was designed based on the nominal model (24). The sampling time for data acquisition was To avoid uncontrollability around θ = 0, the input was decomposed as u =ū + u, and the difference u from the equilibrium inputū was considered as a new input to identify the lifted bilinear model. See Appendix 1 for a discussion of the effect on controllability due to the difference in the input choice.
With the identified lifted bilinear model, the prediction accuracy was evaluated using test data, and the results are shown in Figure 2. For the test data, 30 random initial states that were different from the training data were used and 2[s] of simulations were performed for each trajectory. The predictions are close to the true values for the most of trajectories, showing that the nonlinearity of the target system is well captured by the identified model.

Without state constraints
We evaluated the proposed method using the lifted bilinear model obtained by the identification. First, we show the simulation results for the case without state constraints.
In the simulation, the initial state was set to [−15, 0, 0, 0, 0, 0] T while the target state was set to the origin to control the drone to move horizontally. The simulation time was 10[s], the sampling time was 0.05[s], the prediction horizon length was 30 steps, and the weight matrices in the objective function J 2 were S = 100 · I 6 , Q = I 6 , R = 0.1 · I 2 , R = 0.001 · I 3 . The input constraints were set to T 1 , T 2 ≥ 0. The effect of noise was not taken into account this time.
We used IPOPT with CasADi [15] as the interface for the NMPC solver, and the MATLAB quadprog function for the LMPC solver for proposed method and the Nominal LMPC. We conducted the simulations with MATLAB 2020b running on a Core i5 CPU@1.20GHz, 24.0 GB memory, Windows 10 PC. Figures 3 and 4 show the results of the comparison of state and input with each method, respectively. The dashed lines in the figures indicate constrains. First, for the state shown in Figure 3, all methods successfully stabilized the target system to the origin. However, it can be seen that the behaviour of the Nominal LMPC using the linear approximation model differs significantly from that of the NMPC using the true model, while the proposed method yields behaviour that are almost similar to that of the NMPC. It can also be seen that in the proposed method, the behaviour is more similar to NMPC when there is correction for prediction error than when there is no correction. In particular, after around 2[s], the proposed method almost matches the NMPC. There is a similar trend for the input shown in Figure 4. However, it can be seen that some differences appear for the proposed method compared to NMPC at the time period just after the beginning of the simulation.
Here, the graph of corrective input in the proposed method is shown in Figure 5. From the graph, we can see that the corrective input is relatively large just after the beginning of the simulation, especially until about 1[s], but after that it decreases to small value and  eventually converges to zero. This is considered to be because there were large error in the prior-predicted states at the start of the simulation and also the large prediction error of the time-varying linear model. Then, as the accuracy of the prior-predicted states improve through iterative calculations, the prediction accuracy of the time-varying linear model improves and the corrective input approaches zero because they are no longer needed.
The performance of each method was quantitatively evaluated by the objective function where N is the number of data points obtained in the simulation period. The evaluation results, normalized with the NMPC value as 1, were 1.90 for the Nominal LMPC and 1.01 for the proposed method. These results indicate that the performance of the Nominal LMPC was deteriorated by 90% against the NMPC, whereas the proposed method can achieve control performance comparable to that of the NMPC.

With state constraints
Next, we show the simulation results for the case where |p z | < 0.1 is set as the state constraints. In this case, the initial state was set to [−10, 0, 0, 0, 0, 0, 0] T , and other simulation conditions were the same as in the previous case. Figures 6 and 7 show the comparison results of the state and input by each method, respectively. Even with the state constrains, the behaviour obtained by the proposed method are closer to that of NMPC than Nominal LMPC. Particularly, a major difference is that the constraints on p z are often violated in the Nominal LMPC, while the proposed method satisfies them. This is considered to be achieved by the better prediction accuracy of the proposed method.
The objective function values normalized by the NMPC value were 1.06 for the Nominlal LMPC and 1.00 for the proposed method. These results show that the proposed method outperforms the Nominal LMPC even in the presence of state constraints, and can   achieve control performance comparable to that of the NMPC.

CPU time
Finally, we discuss the results of the computation time evaluation. Figure 8 shows the execution time of the optimization solver for each method for the prediction horizons N p = 10, 20, 30. Figures 8(a) and 8(b) show the evaluation results without and with the state constraints, respectively.
The results show that the proposed method significantly reduces the computation time compared to NMPC, regardless of the presence of the state constraints. On the other hand, the proposed method takes longer computation time than the Nominal LMPC. Main reasons for this are considered to be that the proposed method requires a higher order model by the lifting and the constraints on the lifted states.
From the above results, it was shown through verification that the proposed method has the comparable control performance to NMPC, but at a lower computational cost than NMPC.

Conclusion
In this study, we proposed an LMPC algorithm that can achieve control performance comparable to NMPC with low computational cost, using the identified lifted bilinear model for the input affine system. In the proposed method, the identified lifted bilinear model is used as the prediction model in MPC, and prediction error correction is performed by introducing constraints based on observables and corrective inputs, and then linearizing the nonlinear prediction model and constraints to achieve accurate prediction in LMPC. Comparison with NMPC and Nominal LMPC in numerical simulation shows the superiority of the proposed method in terms of control performance and computational cost.
The proposed method can accelerate time-consum ing optimization calculations in MPC calculations, which is expected to improve the real-time control performance of the actual plants. Although the effectiveness of the proposed method was verified by simulation using a simple model, designing a controller for a more practical plant and verifying it on the real environment remains a future task. In addition, it is also a future challenge to apply the parallel processing method [16] to reduce the computation time of NMPC.
Furthermore, although in this paper the lifted bilinear model was identified assuming the target plant as unknown, in many practical cases a simplified nominal model can be available. So, it is also planned to consider more actively utilizing the nominal model knowledge to extend the proposed method.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Notes on contributors
Masaki Kanai received the B.E. and M.E. degrees from Tokyo Institute of Technology in 2009 and 2011, respectively. He is currently a researcher in the Autonomous Driving Systems Research Department, Controls and Robotics Innovation Center, R&D Group of Hitachi, Ltd., as well as a doctoral student of the Department of Systems and Control Engineering, School of Engineering, Tokyo Institute of Technology. His research interests include nonlinear control, system identification and robotics.
Masaki Yamakita received the B.E., M.E., and Ph.D. degrees from Tokyo Institute of Technology in 1984, 1986, and 1989, respectively. From 1989 he was a research associate in the department of control engineering of the institute, and from 1993 he was a lecturer at Toyohashi University of Technology. He is currently an associate professor in the Department of Systems and Control Engineering of Tokyo Institute of Technology. His research interests include robotics, learning control, and robust and nonlinear control. He is a member of IEEE and RSJ.