Echo State Network Based Model Predictive Control for Active Vibration Control of Hybrid Electric Vehicle Powertrains

Abstract

Reservoir computing refers to a computational framework based on recurrent neural networks that can process time-series data. In an echo state network (ESN), which is a type of reservoir computing framework, the reservoir consists of a recursive network of artificial neurons with nonlinear activation functions. A model predictive control (MPC) technique can determine the control signals by solving the optimization problem of a system using the finite-time domain of each control period. However, real-time optimization cannot be achieved unless the optimal control problem can be solved within the next control period. To overcome this limitation, we propose a new control method based on MPC that explicitly incorporates the predicted disturbance of a time-varying trajectory using ESN to achieve the active vibration control of hybrid electric vehicle (HEV) powertrains. Once the ESN has been trained, the associated MPC explicitly satisfies the constraints over a moving horizon without further training. Instead of completing the real-time optimization within the control period, ESN predicts the future disturbance and applies it to the MPC in the future control period. Based on the predicted future disturbance, the system calculates the optimal control signals required for the future. Thus, real-time control can be realized because the optimal signals are determined before the subsequent control period occurs. The proposed method can be implemented in MPC even if the control period is too short to optimize as long as the disturbance can be reasonably measured and predicted. In this study, the simulation approach was demonstrated using the engine start condition in an HEV powertrain. The importance of this study is that the limitation of MPC relevant to real-time optimization can be relaxed by applying our proposed method.

1. Introduction

Owing to the current climate change crisis, the development of automobiles is encountering several environmental issues. Particularly, it is challenging to satisfy customer demands for driving comfort with superior vehicle performance while ensuring low emissions and high energy efficiency. Moreover, the restrictions pertaining to CO${}_2$ emissions of automobiles have gradually become stricter worldwide. The use of hybrid electric vehicles (HEVs) represents a realistic and promising solution to these problems until the next generation of low emission applications, such as fuel cells, are developed because such vehicles utilize the existing infrastructure. HEVs generate mechanical power from an internal combustion engine (ICE) and electric motors. ICE converts the fuel combustion energy into the driving force. The reciprocating motion of the piston is converted to rotational motion using a crankshaft. However, owing to the periodic combustion and complex components of the drivetrain, ICEs generate torsional oscillations of multiple orders. Although HEVs have improved fuel efficiency and reduced emissions, the ICE is frequently switched on/off while driving. As these operations occur without driver input, the sudden operation may make certain passengers uncomfortable. This discomfort can be primarily attributed to the torsional oscillation of the engine torque. The sudden changes in the driving torque can be alleviated by reducing the oscillation transmitted to the powertrain. Moreover, alternative techniques to reduce fuel consumption and emissions include downsizing and downspeeding. In the downsizing approach, the engine displacement and/or the number of cylinders is reduced, and the associated performance deterioration is compensated via turbocharging. Although energy loss can be reduced in this manner, the vibration may intensify owing to the high combustion pressure and low oscillation speed, which is undesirable.

To achieve active vibration reduction of HEV powertrains by focusing on torque oscillation during engine start, this study proposes a motor control method using model predictive control (MPC). MPC is a control method that can determine the control input by solving the optimal control problem for a finite future at every control period. Because optimal control can be realized while satisfying the constraints, MPC is advantageous in various applications. However, it has been observed that real-time optimization cannot be realized unless the optimal control problem is solved numerically within the sampling period [1]. In the proposed approach, to facilitate the realization of real-time optimization, we determined the control input for the MPC to be applied in the future state in advance instead of conducting the optimization within the next control period. Specifically, we attempted to improve the control performance of MPC by using an echo state network (ESN) to learn the time trajectory of past disturbances. Additionally, by predicting future disturbances pertaining to the MPC, we aim to achieve active vibration control of the HEV powertrain.

The remaining paper is organized as follows: Section 2 describes the research background and state-of-the-art solution techniques. Section 3 and Section 4 describe the proposed control approach and simulation results for the case study, respectively. Section 5 discusses the results for the considered implementation, and Section 6 presents the concluding remarks.

2. Related Work

ICEs induce torsional vibrations, which may cause passenger discomfort. Such torsional vibrations can be reduced via passive and active methods. As a passive approach, spring-mass dampers are widely used in existing ICE vehicles. Effective damping can be achieved by tuning the spring rate and inertia of the components. However, downsizing and downspeeding can cause low-frequency vibrations with a large amplitude during ICE operation at speeds slightly higher than the idle speed. Moreover, typical conventional dampers do not exhibit a satisfactory damping performance under such operating conditions. Since 2008, the use of centrifugal pendulum-type absorbers (CPAs) has increased as an alternative passive approach [2]. A CPA isolates the torsional vibrations between the ICE and transmission. According to the physical principle of pendulum system dynamics, the order of excitation and pendulum vibrations are balanced. Hence, the target order of the vibration can be compensated regardless of the frequency range related to the engine speed. In contrast to passive approaches, active approaches can satisfy the low-cost and high-performance requirements by maximizing the capability of motor control in HEVs. Numerous approaches have been proposed to realize active vibration reduction of powertrains. For example, a reversible alternator, as an active flywheel, was developed to replace the flywheel generally used to reduce engine speed oscillation [3]. Additionally, a harmonic activation neural network was applied to reduce vibrations by using a starter generator to apply an inverse torque oscillation to the crankshaft to enable adaptive damping [4]. Furthermore, as a practical approach, explicit dead-time compensation was implemented in a control scheme to address the lag behavior of the actuator and time delay [5]. Pham et al. [6] proposed an anti-jerk controller to reduce the undesired drivetrain oscillations for an HEV with a feedforward and feedback strategy involving dead-time compensation. In this approach, the sampling time of the control period was set as 10 ms, whereas the dead-time delay was set to be approximately twice the control period. Although the contribution of dead-time compensation in vibration reduction control has been demonstrated, the control performance in cases where the dead-time is several times larger than the sampling period is unclear. Moreover, Vadamalu et al. [7] applied a time-delayed input in the MPC and demonstrated time-variant constraint handling by using the explicit dead-time compensation scheme to achieve active vibration reduction of HEV powertrains. In this approach, the unknown future disturbance up to the predictive horizon was predicted by accumulating the system matrix of the state-space representation by utilizing the explicit dead-time scheme. However, the presence of future disturbances, which often vary significantly with time, in real systems was not considered. The existing MPC approaches employing the explicit dead-time may lack precision in such cases because the error of the disturbance between the predicted and actual values is larger than that in the case of stationary or cyclic disturbance conditions.

To create MPC prediction models, neural networks and several other machine learning techniques have been successfully applied in various types of control applications [8,9,10]. Pan et al. [11] adopted recurrent neural networks with input-output training to realize the unknown nonlinear system identification. Furthermore, they decomposed the non-convex optimization problem associated with a nonlinear MPC via Taylor expansion at each time step. Some researchers replaced existing approaches involving recurrent neural networks wiith ESNs to achieve predictive control with improved learning efficiency. Furthermore, ESNs have been applied successfully in the MPC framework. Xiang et al. [12] applied the Taylor expansion to linearize the ESN in the context of an operating point and adopted an ESN to compensate for the truncation error in a partially observed dynamical system. Jordanou et al. [13] implemented ESNs into the practical nonlinear MPC framework developed by Plucenio et al. [14] by applying the analytically computed gradient from the ESN model to the predictive model. Zhang et al. [15] applied an ESN for the decentralized control problem of continuous-time nonlinear interconnected systems. They treated the interconnected terms between the subsystems as the disturbances added to the system dynamics. Liu et al. [16] presented an improved algorithm to solve the optimal control problem of nonlinear systems with the Hamilton–Jacobi–Bellman (HJB) equation approximated by an ESN to solve the optimal control problem. In inference-based modeling, explicit knowledge regarding the structure of the mathematical model is not required. Moreover, highly nonlinear dynamics under uncertain disturbances can be modeled. Nevertheless, although inference-based modeling is useful to realize system identification and control, a clear relationship does not exist between the target system parameters and the associated effects, such as grey-box modeling. Furthermore, in real-time control, considerable time is required to realize MPC calculation within a control period. Considering the aforementioned aspects, to overcome the limitations of inference-based models and to enhance the model performance under uncertain disturbances, we propose an ESN-based MPC to achieve active vibration control of HEV powertrains. The key contributions of this study are as follows:

• A novel control scheme based on an ESN and MPC is developed to improve the control performance in cases involving highly irregular disturbances.
• The proposed MPC control scheme can be applied even if the optimization calculation does not finish within the control period.
• The predictive model for the HEV system is based on the state-space representation to understand the relationship between the model parameters and output.
• The proposed method can be applied to several other control targets apart from the active vibration control of HEV powertrains.

3. Control Method

This section describes the theory of the developed ESN-based MPC scheme. The control target was an HEV powertrain with two motors and one planetary gear set, as illustrated in Figure 1. The ICE and two motors (M1 and M2) were connected with the planetary gear set. The driving force generated by the ICE and motors was distributed to the output shaft of the drivetrain. The torque at the engine start was simulated according to the engine simulator, which was provided as a sample model [17]. The output of the simulator included the in-cylinder pressure of each cylinder and rotation angle of the crankshaft. The output torque of the crankshaft was calculated using a slider-crank mechanism.

This section is structured as follows. Section 3.1 presents an overview of the proposed method. Section 3.2 describes the basic equations of the target HEV system. The theoretical treatment to solve the constrained problem using the MPC is described in Section 3.3 and Section 3.4. To compare the proposed approach with a conventional approach, an overview of the explicit dead-time compensation scheme is presented in Section 3.5. The method to predict the time-varying disturbance by using the ESN is described in Section 3.6.

3.1. Concept of ESN-MPC

In real-life circumstances of the control target, the system state changes dynamically, and unexpected events may occur. MPC, as a control method, can help predict future events and optimize control signals at each control period. The objective of the MPC is to control and optimize system performance for a predicted finite-time future. However, to realize real-time optimization, the MPC generally requires full calculation results within the control period. Therefore, this framework cannot be easily applied in a system where the control period is too short to complete the calculation. Instead of performing real-time optimization within the subsequent control period, the forecasted future disturbance is used to establish the future optimal control. The ESN predicts the future disturbance, which is multiple steps ahead of the period, to calculate the optimal control signals required by the MPC in the future. Figure 2 shows the controller structure in the proposed method. The upper left part illustrates the prediction sequence of a future disturbance in the system. The difference in time between the current calculation and optimization calculation is termed the time margin. The result predicted using the ESN is applied to the MPC in the future to achieve optimization at the time when the control signals are required. In other words, predicting the future disturbance provides a time margin to achieve the optimization calculation.

3.2. Basic Equations for the HEV Powertrain

The control target, i.e., the HEV drivetrain system, was configured to have a torsional damper, a planetary gear set, and two electric motors. The engine torque was input as a system disturbance. The equation of the system’s motion was derived using the Lagrangian equation. Considering the characteristics of the planetary gears, the relational expression for the rotational speed can be defined as follows:

$$\begin{array}{c}\hfill {\omega }_{m2}=(\lambda +1){\omega }_c-\lambda {\omega }_{m1},{\omega }_p=-\frac{\lambda +1}{\lambda -1}{\omega }_c+\frac{2\lambda }{\lambda -1}{\omega }_{m1}\end{array}$$

(1)

where $\lambda$ and $\omega$ denote the planetary gear ratio and rotation speed, respectively. The subscripts c, $m1$, $m2$, and p correspond to the carrier, motor 1, motor 2, and pinion gear, respectively. The relationship between each torque can be expressed as follows:

$$\begin{array}{c}\hfill T_r=\frac{T_{m1}}{\lambda }+T_{m2},T_c=(\lambda +1)T_r.\end{array}$$

(2)

where $T_r$, $T_c$, $T_{m1}$, and $T_{m2}$ indicate the torque of the ring gear, carrier, and motors 1 and 2, respectively. The kinetic energy of the entire system $E_k$, potential energy owing to the damper mechanism $E_u$, and dissipation energy $E_d$ can be expressed as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill E_k=\frac{1}{2}\{I_e{\omega }_e^2+I_c{\omega }_c^2+I_{m1}{\omega }_{m1}^2+I_{m2}{\omega }_{m2}^2+I_p{\omega }_p^2+n_p{m}_p{(R_s+R_p)}^2{\omega }_c^2\},\end{array}{E}_u=\frac{1}{2}{k}_{dp}{({\theta }_e-{\theta }_c)}^2,E_d=\frac{1}{2}{c}_{dp}{({\omega }_e-{\omega }_c)}^2\end{array}$$

(3)

where $\theta$ is the rotation angle; I is the moment of inertia; $k_{dp}$ and $c_{dp}$ denote the torsional stiffness and damping coefficient of the damper, respectively; $R_s$ and $R_p$ denote the sun gear and pinion gear radius, respectively; $m_p$ is the pinion gear mass; and $n_p$ denotes the number of pinion gears. The subscript e used with $\theta$, I, and $\omega$ indicates the engine crankshaft. In the Lagrangian equation, the Lagrangian is expressed as $L=E_k$ − $E_u$, the dissipation energy is denoted as $D=E_d$, and the external force in this system is represented as $f_i$. The rotational speed of each element is as follows: ${\omega }_e={\dot{\theta }}_e$, ${\omega }_c={\dot{\theta }}_c$, ${\omega }_{m1}={\dot{\theta }}_{m1}$, ${\omega }_{m2}={\dot{\theta }}_{m2}$, and ${\omega }_p={\dot{\theta }}_p$. The Lagrangian equation can be expressed as follows:

$$\begin{array}{c}\hfill \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{{\theta }_i}}\right)-\frac{\partial L}{\partial {\theta }_i}+\frac{\partial D}{\partial {\dot{\theta }}_i}=f_i(i=1,2,3).\end{array}$$

(4)

In the planetary gear system, ${\omega }_{m1}$ and ${\omega }_p$ were eliminated using the relationship from Equation (1). The basic equations were implemented as a linear system with three degrees of freedom. The rotational speed and angle of each element i in Equation (4) are set as follows: ${\dot{\theta }}_1={\omega }_e$, ${\dot{\theta }}_2={\omega }_c$, ${\dot{\theta }}_3={\omega }_{m2}$, ${\theta }_1={\theta }_e$, ${\theta }_2={\theta }_c$, and ${\theta }_3={\theta }_{m2}$. By solving Equations (1)–(4) for $i=1$–3, we obtain Equation (5), which is the equation of the system’s motion:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{I}_e{\dot{\omega }}_e\\ J_1{\dot{\omega }}_c\\ J_1{\dot{\omega }}_{m2}\end{array}\right]=\left[\begin{array}{cc}{k}_{dp}& c_{dp}\\ J_2{k}_{dp}& J_2{c}_{dp}\\ J_3{k}_{dp}& J_3{c}_{dp}\end{array}\right]\left[\begin{array}{c}{\theta }_e-{\theta }_c\\ {\omega }_e-{\omega }_c\end{array}\right]+\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\end{array}\right]\end{array}$$

(5)

where

$$\begin{array}{cc}& f_1=T_e,f_2=2{(\lambda +1)}^2{J}_2\left(\frac{T_{m1}}{\lambda }+T_{m2}\right),f_3=2{(\lambda +1)}^2{J}_4\left(\frac{T_{m1}}{\lambda }+T_{m2}\right),\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{cc}& J_1=J_2{(R_s+R_p)}^2{n}_p{m}_p+\left\{{({\lambda }^2-2\lambda -1)}^2{I}_c+{(\lambda +1)}^2(I_r+{(\lambda -2)}^2{I}_s)\right\}{I}_p{n}_p+{(\lambda -1)}^2\left\{I_c(I_s+{\lambda }^2{I}_r)+{(\lambda +1)}^2{I}_r{I}_s\right\},\hfill \end{array}\end{array}$$

$$\begin{array}{cc}& J_2={({\lambda }^2-2\lambda -1)}^2{I}_p{n}_p+{(\lambda -1)}^2(I_s+{\lambda }^2{I}_r),J_3=-(\lambda +1)J_4,J_4=({\lambda }^2-2\lambda -1)I_p{n}_p-{(\lambda -1)}^2{I}_s.\end{array}$$

Because the torque transmitted by the torsional damper is generated by the relative revolution of the engine crankshaft and carrier, the angle and angular velocity is presented as ${\theta }_e$ − ${\theta }_c$ and ${\omega }_e$ − ${\omega }_c$, according to Equation (5).

Considering a dynamical system with some delay in the actuation, our control scheme provides future control signals in which future disturbances are predicted using ESN.

3.3. State-Space Representation and Kalman Filter

The target system is expressed as a linear time-invariant system with the state-space representation, as follows:

$$\begin{array}{ccc}& & \dot{\mathit{x}}={\mathit{A}}_c\mathit{x}+{\mathit{B}}_c\mathit{u}+{\mathit{B}}_{ce}\mathit{w},\hfill \end{array}$$

(6)

$$\begin{array}{ccc}& & \mathit{y}={\mathit{C}}_c\mathit{x}\hfill \end{array}$$

(7)

with

$$\begin{array}{ccc}& & \mathit{x}={[{\theta }_e-{\theta }_c,{\omega }_e-{\omega }_c,{\omega }_{m2}]}^{\mathrm{T}},\mathit{y}={[{\omega }_{m2},T_r]}^{\mathrm{T}},\mathit{u}={[T_{m1},T_{m2}]}^{\mathrm{T}},\mathit{w}=T_e,\hfill \end{array}$$

$$\begin{array}{ccc}& & {\mathit{A}}_c=\left[\begin{array}{ccc}0& 1& 0\\ \left(\frac{1}{I_e}-\frac{J_2}{J_1}\right)k_{dp}& \left(\frac{1}{I_e}-\frac{J_2}{J_1}\right)c_{dp}& 0\\ \frac{J_3}{J_1}{k}_{dp}& \frac{J_3}{J_1}{c}_{dp}& 0\end{array}\right],{\mathit{B}}_c=\frac{2{(\lambda +1)}^2}{\lambda }\left[\begin{array}{cc}0& 0\\ J_2& \lambda J_2\\ J_4& \lambda J_4\end{array}\right],{\mathit{B}}_{ce}={[0,1,0]}^{\mathrm{T}},{\mathit{C}}_c=\left[\begin{array}{ccc}0& 0& 1\\ \frac{k_{dp}}{\lambda +1}& \frac{c_{dp}}{\lambda +1}& 0\end{array}\right]\hfill \end{array}$$

where $\mathit{x}$, $\mathit{y}$, $\mathit{u}$, and $\mathit{w}$ are the system state, system output, control input, and disturbance, respectively. The engine torque $T_e$ is treated as a disturbance acting on the system, and the control inputs, $T_{m1}$ and $T_{m2}$, are subject to the following constraints and limit values:

$$\begin{array}{ccc}& & T_{m1}^{min}\le T_{m1}\le T_{m1}^{max},T_{m2}^{min}\le T_{m2}\le T_{m2}^{max}.\hfill \end{array}$$

(8)

Furthermore, it was assumed that the target system contains elements that cannot be observed directly, and a discrete time-varying Kalman filter [18] was embedded to estimate the state with noise. The equations of the continuous time system (Equations (6) and (7)) can be discretized and modified by adding process noise ${\mathit{n}}_k$ and output noise ${\mathit{v}}_k$, as follows:

$$\begin{array}{ccc}\hfill {\mathit{x}}_{k+1}& =& {\mathit{A}\mathit{x}}_k+{\mathit{B}\mathit{u}}_k+{\mathit{B}}_e{\mathit{w}}_k+{\mathit{n}}_k,\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\mathit{y}}_k& =& {\mathit{C}\mathit{x}}_k+{\mathit{v}}_k\hfill \end{array}$$

(10)

where $\mathit{A}$, $\mathit{B}$, ${\mathit{B}}_e$, and $\mathit{C}$ are obtained from the discrete time state space representation. k represents a time step in discrete time. The algorithm of discrete time-varying Kalman filter is given by the following equations:

$$\begin{array}{ccc}\hfill {\widehat{\mathit{x}}}_{k|k}& =& {\widehat{\mathit{x}}}_{k|k-1}+{\mathit{K}}_k({\mathit{y}}_k-\mathit{C}{\widehat{\mathit{x}}}_{k|k-1}),\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill {\widehat{\mathit{x}}}_{k|k-1}& =& \mathit{A}{\widehat{\mathit{x}}}_{k-1|k-1}+{\mathit{B}\mathit{u}}_{k-1}+{\mathit{B}}_e{\mathit{w}}_{k-1},\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill {\widehat{\mathit{y}}}_{k|k-1}& =& \mathit{C}{\widehat{\mathit{x}}}_{k|k-1}\hfill \end{array}$$

(13)

where ${\widehat{\mathit{x}}}_k:={\widehat{\mathit{x}}}_{k|k}$ denotes the estimated state at the time step k. Now, the covariance matrices are set to ${\mathit{Q}}_k$ and ${\mathit{R}}_k$, regarding each of ${\mathit{n}}_k$ and ${\mathit{v}}_k$. Then, the Kalman gain ${\mathit{K}}_k$ can be updated in the following manner:

$$\begin{array}{ccc}\hfill {\mathit{K}}_k& =& {\mathit{P}}_{k|k-1}{\mathit{C}}^{\mathrm{T}}{{\mathit{S}}_k}^{-1},\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill {\mathit{P}}_{k|k-1}& =& \mathit{A}{\mathit{P}}_{k-1|k-1}{\mathit{A}}^{\mathrm{T}}+{\mathit{Q}}_k,\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill {\mathit{S}}_k& =& \mathit{C}{\mathit{P}}_{k|k-1}{\mathit{C}}^{\mathrm{T}}+{\mathit{R}}_k,\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill {\mathit{P}}_{k|k}& =& (\mathit{I}-{\mathit{K}}_k\mathit{C}){\mathit{P}}_{k|k-1}.\hfill \end{array}$$

(17)

3.4. Solving the Predictive Control Problem

The predicted state for next time step ${\widehat{\mathit{x}}}_{k+1}$ was defined based on the current estimated state ${\widehat{\mathit{x}}}_k$ filtered by Kalman filter and using the given predicted values ${\widehat{\mathit{u}}}_k$ and ${\widehat{\mathit{w}}}_k$ as follows:

$$\begin{array}{c}\hfill {\widehat{\mathit{x}}}_{k+1}={\mathit{A}\widehat{\mathit{x}}}_k+{\mathit{B}\widehat{\mathit{u}}}_k+{\mathit{B}}_e{\widehat{\mathit{w}}}_k,{\widehat{\mathit{y}}}_k={\mathit{C}\widehat{\mathit{x}}}_k.\end{array}$$

(18)

The predicted control input ${\widehat{\mathit{u}}}_k$ can be expressed as the sum of the predicted input difference $\mathsf{\Delta }{\widehat{\mathit{u}}}_k$ and control input at the previous step ${\mathit{u}}_{k-1}$ as follows:

$$\begin{array}{ccc}\hfill {\widehat{\mathit{u}}}_k& =& \mathsf{\Delta }{\widehat{\mathit{u}}}_k+{\mathit{u}}_{k-1}.\hfill \end{array}$$

(19)

${\widehat{\mathit{w}}}_k$ can be estimated using the ESN in the proposed ESN-MPC. Conversely, a conventional EDT-MPC follows the method that is described in Section 3.5. Assume that the variable $\mathsf{\Delta }{\mathit{U}}_k={[\mathsf{\Delta }{\widehat{\mathit{u}}}_k,\cdots ,\mathsf{\Delta }{\widehat{\mathit{u}}}_{k+H_u-1}]}^{\mathrm{T}}$ represents the optimal control sequence at the time step k. $\mathsf{\Delta }{\widehat{\mathit{u}}}_k$ in Equation (19) is derived by Equation (22) as the control sequence $\mathsf{\Delta }{\mathit{U}}_k$, which is described later. The predicted output ${\mathit{Y}}_k$ can be expressed as follows [19]:

$$\begin{array}{c}\hfill {\mathit{Y}}_k={\mathsf{\Psi }\widehat{\mathit{x}}}_k+{\mathit{\gamma }\mathit{u}}_{k-1}+\mathsf{\Theta }\mathsf{\Delta }{\mathit{U}}_k+{\mathsf{\Xi }{\mathit{W}}_{\mathit{m}}}_k\end{array}$$

(20)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathit{Y}}_k=\left[\begin{array}{c}{\widehat{\mathit{y}}}_{k+1}\\ ⋮\\ {\widehat{\mathit{y}}}_{k+H_u}\\ {\widehat{\mathit{y}}}_{k+H_u+1}\\ ⋮\\ {\widehat{\mathit{y}}}_{k+H_p}\end{array}\right]& \hfill ,{{\mathit{W}}_{\mathit{m}}}_k=\left[\begin{array}{c}{\widehat{\mathit{w}}}_k\\ ⋮\\ {\widehat{\mathit{w}}}_{k+H_u-1}\\ {\widehat{\mathit{w}}}_{k+H_u}\\ ⋮\\ {\widehat{\mathit{w}}}_{k+H_p-1}\end{array}\right],\end{array}\begin{array}{cc}\hfill \mathsf{\Psi }=\left[\begin{array}{c}\mathit{C}\mathit{A}\\ ⋮\\ {\mathit{C}\mathit{A}}^{H_u}\\ {\mathit{C}\mathit{A}}^{H_u+1}\\ ⋮\\ {\mathit{C}\mathit{A}}^{H_p}\end{array}\right]& \hfill ,\mathit{\gamma }=\left[\begin{array}{c}\mathit{C}\mathit{B}\\ ⋮\\ \sum _{i=0}^{H_u-1}{\mathit{C}\mathit{A}}^i\mathit{B}\\ \sum _{i=0}^{H_u}{\mathit{C}\mathit{A}}^i\mathit{B}\\ ⋮\\ \sum _{i=0}^{H_p-1}{\mathit{C}\mathit{A}}^i\mathit{B}\end{array}\right],\end{array}\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathsf{\Theta }=\left[\begin{array}{cccc}\mathit{C}\mathit{B}& 0& \cdots & 0\\ \mathit{C}\mathit{A}\mathit{B}+\mathit{C}\mathit{B}& \mathit{C}\mathit{B}& \cdots & 0\\ ⋮& \cdots & \ddots & ⋮\\ \sum _{i=0}^{H_u-1}{\mathit{C}\mathit{A}}^i\mathit{B}& \cdots & \cdots & \mathit{C}\mathit{B}\\ \sum _{i=0}^{H_u}{\mathit{C}\mathit{A}}^i\mathit{B}& \cdots & \cdots & \mathit{C}\mathit{A}\mathit{B}+\mathit{C}\mathit{B}\\ ⋮& \cdots & \ddots & ⋮\\ \sum _{i=0}^{H_p-1}{\mathit{C}\mathit{A}}^i\mathit{B}& \sum _{i=0}^{H_p-2}{\mathit{C}\mathit{A}}^i\mathit{B}& \cdots & \sum _{i=0}^{H_p-H_u}{\mathit{C}\mathit{A}}^i\mathit{B}\end{array}\right],\end{array}\mathsf{\Xi }=\left[\begin{array}{cccc}{\mathit{C}\mathit{B}}_e& \cdots 0& \cdots & 0\\ {\mathit{C}\mathit{A}\mathit{B}}_e& {\mathit{C}\mathit{B}}_e& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ {\mathit{C}\mathit{A}}^{H_p-1}{\mathit{B}}_e& {\mathit{C}\mathit{A}}^{H_p-2}{\mathit{B}}_e& \cdots & {\mathit{C}\mathit{B}}_e\end{array}\right]\end{array}$$

where $H_u$ and $H_p$ represent the control horizon and prediction horizon, respectively, and $H_u\le H_p$. The MPC, which is based on the constrained finite-horizon optimization, for the tracking problem can be formulated as follows:

$$\begin{array}{c}\hfill min{J}_k=\sum _{i=0}^{H_p}\Vert {\widehat{\mathit{y}}}_{k+i}-{{\mathit{r}}_{\mathit{t}}}_{k+i}{\Vert }_{{}_{\mathcal{Q}}}^2+\sum _{i=0}^{H_u}{\Vert \mathsf{\Delta }{\widehat{\mathit{u}}}_{k+i}\Vert }_{{}_{\mathcal{R}}}^2\end{array}$$

(21)

subject to

$$\begin{array}{ccc}& & {\widehat{\mathit{x}}}_{k+i+1}={\mathit{A}\widehat{\mathit{x}}}_{k+i}+{\mathit{B}\widehat{\mathit{u}}}_{k+i}+{\mathit{B}}_e{\widehat{\mathit{w}}}_{k+i},{\widehat{\mathit{y}}}_{k+i}={\mathit{C}\widehat{\mathit{x}}}_{k+i},{\widehat{\mathit{x}}}_0={\mathit{x}}_0,\hfill \\ & & {\widehat{\mathit{y}}}^{min}\le {\widehat{\mathit{y}}}_{k+i}\le {\widehat{\mathit{y}}}^{max},{\widehat{\mathit{u}}}^{min}\le {\widehat{\mathit{u}}}_{k+i}\le {\widehat{\mathit{u}}}^{max},\mathsf{\Delta }{\widehat{\mathit{u}}}_{k+i}=0fori\in \{H_u+1,\cdots ,H_p\}\hfill \end{array}$$

where the notations ${\Vert ·\Vert }_{{}_{\mathcal{Q}}}^2$ and ${\Vert ·\Vert }_{{}_{\mathcal{R}}}^2$ represent the weighted norms as a quadratic form, which indicates, e.g., for $\mathcal{Q}$, that ${\mathit{x}}^{\mathrm{T}}\mathcal{Q}\mathit{x}=\Vert {\mathcal{Q}}^{1/2}{\mathit{x}\Vert }^2={\Vert \mathit{x}\Vert }_{{}_{\mathcal{Q}}}^2$. Here, i is the step number, $\widehat{\mathit{y}}$ is the predicted output, ${\mathit{r}}_{\mathit{t}}$ is the reference trajectory for the output, and $\widehat{\mathit{w}}$ is the predicted disturbance. The diagonal matrices $\mathcal{Q}$ and $\mathcal{R}$ represent the output and control input weights, respectively. The standard predictive control problem can be solved by updating the solution each time. This modification is performed each time with optimal solution changes. This type of method is called real-time iteration [20]. Herein, we consider the target HEV powertrain as a linear time-invariant system. The predictive control problem of this system involves certain constraints. To solve the constrained optimization problem, it can be formulated as a quadratic programming (QP) problem based on the discretization method in Equations (18)−(21), using a cost function and predicted output [19].

$$\begin{array}{c}\hfill \underset{\mathsf{\Delta }{U}_k}{min}\frac{1}{2}\mathsf{\Delta }{\mathit{U}}_k^{\mathrm{T}}\mathsf{\Phi }\mathsf{\Delta }{\mathit{U}}_k-{\mathcal{G}}^{\mathrm{T}}\mathsf{\Delta }{\mathit{U}}_k\end{array}$$

(22)

subject to

$$\begin{array}{c}\hfill \mathsf{\Omega }\mathsf{\Delta }{\mathit{U}}_k\le \mathit{\eta }\end{array}$$

where

$$\begin{array}{c}\hfill \mathcal{G}=2{\mathsf{\Theta }}^{\mathrm{T}}\mathcal{Q}{\mathit{\epsilon }}_k,\mathsf{\Phi }={\mathsf{\Theta }}^{\mathrm{T}}\mathcal{Q}\mathsf{\Theta }+\mathcal{R},\end{array}$$

$$\begin{array}{c}\hfill \mathsf{\Omega }=\left[\begin{array}{c}-{\mathit{C}}_{\mathbf{2}}\\ {\mathit{C}}_{\mathbf{2}}\\ -\mathsf{\Theta }\\ \mathsf{\Theta }\end{array}\right],\mathit{\eta }=\left[\begin{array}{c}-{\mathit{u}}^{min}+{\mathit{C}}_{\mathbf{1}}{\mathit{u}}_{k-1}\\ {\mathit{u}}^{max}-{\mathit{C}}_{\mathbf{1}}{\mathit{u}}_{k-1}\\ -{\mathit{y}}^{min}+{\mathsf{\Psi }\mathit{x}}_k+{\mathit{\gamma }\mathit{u}}_{k-1}+{\mathsf{\Xi }{\mathit{W}}_{\mathit{m}}}_k\\ {\mathit{y}}^{max}-{\mathsf{\Psi }\mathit{x}}_k-{\mathit{\gamma }\mathit{u}}_{k-1}-{\mathsf{\Xi }{\mathit{W}}_{\mathit{m}}}_k\end{array}\right],{\mathit{C}}_{\mathbf{1}}=\left[\begin{array}{c}\mathit{I}\\ \mathit{I}\\ \mathit{I}\\ ⋮\\ \mathit{I}\end{array}\right],{\mathit{C}}_{\mathbf{2}}=\left[\begin{array}{ccccc}\mathit{I}& 0& 0& \cdots & 0\\ \mathit{I}& \mathit{I}& 0& \cdots & 0\\ \mathit{I}& \mathit{I}& \mathit{I}& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ \mathit{I}& \mathit{I}& \mathit{I}& \cdots & \mathit{I}\end{array}\right],\end{array}$$

$$\begin{array}{cc}& {\mathit{\epsilon }}_k={\mathcal{T}}_k-{\mathsf{\Psi }\mathit{x}}_k-{\mathit{\gamma }\mathit{u}}_{k-1}-{\mathsf{\Xi }{\mathit{W}}_{\mathit{m}}}_k,{\mathcal{T}}_k=\left[\begin{array}{c}{{\mathit{r}}_{\mathit{t}}}_{k+1|k}\\ ⋮\\ {{\mathit{r}}_{\mathit{t}}}_{k+H_p|k}\end{array}\right]\end{array}$$

where ${\mathit{\epsilon }}_k$ is the tracking error that changes in each step owing to the system disturbance. ${\mathcal{T}}_k$ is the reference matrix. ${\mathit{u}}^{min}$ and ${\mathit{u}}^{max}$ denote the input limits, and ${\mathit{y}}^{min}$ and ${\mathit{y}}^{max}$ denote the output limits of the constraints. Owing to the generalized formulation, standard algorithms for the QP problem can be applied to solve the considered problem.

3.5. Explicit Dead-Time Compensation-Based MPC (EDT-MPC)

We compared the proposed and a conventional method in MPC. The explicit dead-time compensation-based MPC (EDT-MPC) was selected as the conventional method. Generally, a controlled system involves a certain dead-time. Additionally, the time lag behavior of the actuator must also be considered. To avoid the presence of high order matrices in the augmented representation, the EDT compensation was developed by Santos et al. [5], who examined this aspect in terms of stability and constraint satisfaction for discrete-time linear systems. The state prediction by EDT represents the current time, including the dead-time. In this study, the state prediction by the EDT was extended to the future state prediction for comparison with the proposed method using ESN.

Consider the following uncertain discrete-time linear system with dead-time:

$$\begin{array}{c}\hfill {\mathit{x}}_{k+1}={\mathit{A}\mathit{x}}_k+{\mathit{B}\mathit{u}}_{k-d}+{\mathit{B}}_e{\mathit{w}}_k\end{array}$$

(23)

where d represents the nominal process dead-time, and ${\mathit{w}}_k$ is the vector of disturbance. In this study, d is treated as an artificial dead-time in EDT-MPC for comparison. The proposed ESN-MPC uses the time margin, which is the time difference between the current and actual execution times (Figure 2). Consider the following prediction model, where the system to be controlled depends only on past control actions:

$$\begin{array}{c}\hfill {\mathit{x}}_{k+1}={\mathit{A}\mathit{x}}_k+{\mathit{B}\mathit{u}}_{k-d}.\end{array}$$

(24)

The predicted value for ${\mathit{x}}_{k+d}$, based on the information available at the k step, denoted by ${\mathit{x}}_{k+d|k}$, can be obtained using Equation (24) recursively, as follows:

$$\begin{array}{ccc}& & {\mathit{x}}_{k+d|k}={\mathit{A}}^d{\mathit{x}}_k+\sum _{j=1}^d\left[{\mathit{A}}^{j-1}{\mathit{B}\mathit{u}}_{k-j}\right],\hfill \end{array}$$

(25)

$$\begin{array}{ccc}& & {\tilde{\mathit{x}}}_k:={\mathit{x}}_{k+d|k}.\hfill \end{array}$$

(26)

By considering the explicit compensation scheme (Santos et al. [5]), given by Equations (25) and (26), the predicted behavior can be noted as

$$\begin{array}{ccc}\hfill {\tilde{\mathit{x}}}_{k+1}& =& \mathit{A}{\tilde{\mathit{x}}}_k+{\mathit{B}\mathit{u}}_k+{\mathit{B}}_e{\tilde{\mathit{w}}}_k\hfill \end{array}$$

(27)

where ${\tilde{\mathit{w}}}_k:={\mathit{w}}_{k+d|k}$ is the effect of ${\mathit{w}}_k$ on the predicted state ${\tilde{\mathit{x}}}_k$, which is uncertain at the timing of the k step. From Equation (27), ${\mathit{B}}_e{\tilde{\mathit{w}}}_k$ can be obtained as:

$$\begin{array}{ccc}\hfill {\mathit{B}}_e{\tilde{\mathit{w}}}_k& =& {\tilde{\mathit{x}}}_{k+1}-\mathit{A}{\tilde{\mathit{x}}}_k-{\mathit{B}\mathit{u}}_k.\hfill \end{array}$$

(28)

By applying Equations (23), (25) and (26), Equation (28) transforms into the following equation:

$$\begin{array}{ccc}\hfill {\mathit{B}}_e{\tilde{\mathit{w}}}_k& =& {\mathit{A}}^d{\mathit{x}}_{k+1}+\sum _{j=1}^d\left[{\mathit{A}}^{j-1}{\mathit{B}\mathit{u}}_{k-j+1}\right]-\mathit{A}\left\{{\mathit{A}}^d{\mathit{x}}_k+\sum _{j=1}^d\left[{\mathit{A}}^{j-1}\mathit{B}{\mathit{u}}_{k-j}\right]\right\}-{\mathit{B}\mathit{u}}_k\hfill \\ & =& {\mathit{A}}^d[{\mathit{x}}_{k+1}-\mathit{A}{\mathit{x}}_k-{\mathit{B}\mathit{u}}_{k-d}]\hfill \\ & =& {\mathit{A}}^d{\mathit{B}}_e{\mathit{w}}_k.\hfill \end{array}$$

(29)

Then, the uncertain vector of ${\mathit{B}}_e{\tilde{\mathit{w}}}_k$ can be determined by ${\mathit{B}}_e{\mathit{w}}_k$. Using Equation (23) recursively and applying ${\mathit{w}}_{k-j+d|k}={\tilde{\mathit{w}}}_{k-j}$, as shown in the following equation:

$$\begin{array}{ccc}\hfill {\mathit{x}}_{k+d|k}& =& {\mathit{A}}^d{\mathit{x}}_k+\sum _{j=1}^d\left[{\mathit{A}}^{j-1}{\mathit{B}\mathit{u}}_{k-j}\right]+\sum _{j=1}^d\left[{\mathit{A}}^{j-1}{\mathit{B}}_e{\mathit{w}}_{k-j+d|k}\right]\hfill \\ & =& {\mathit{A}}^d{\mathit{x}}_k+\sum _{j=1}^d\left[{\mathit{A}}^{j-1}{\mathit{B}\mathit{u}}_{k-j}\right]+\sum _{j=1}^d\left[{\mathit{A}}^{j-1}{\mathit{B}}_e{\tilde{\mathit{w}}}_{k-j}\right].\hfill \end{array}$$

(30)

By incorporating Equation (29) into Equation (30), the predicted value for ${\mathit{x}}_{k+d|k}$ can be obtained as

$$\begin{array}{c}\hfill {\mathit{x}}_{k+d|k}={\mathit{A}}^d{\mathit{x}}_k+\sum _{j=1}^d\left[{\mathit{A}}^{j-1}{\mathit{B}\mathit{u}}_{k-j}\right]+\sum _{j=1}^d\left[{\mathit{A}}^{d+j-1}{\mathit{B}}_e{\mathit{w}}_{k-j}\right].\end{array}$$

(31)

This extended representation does not include any unknown information available at the k step. To implement the conventional EDT-MPC, the problem formulation in Equation (21) can be extended as follows:

$$\begin{array}{c}\hfill min{J}_k=\sum _{i=0}^{H_p}\Vert {\widehat{\mathit{y}}}_{k+d+i|k}-{{\mathit{r}}_{\mathit{t}}}_{k+d+i}{\Vert }_{{}_{\mathcal{Q}}}^2+\sum _{i=0}^{H_u}{\Vert \mathsf{\Delta }{\widehat{\mathit{u}}}_{k+i}\Vert }_{{}_{\mathcal{R}}}^2\end{array}$$

(32)

subject to

$$\begin{array}{ccc}& & {\widehat{\mathit{x}}}_{k+d+i+1|k}={\mathit{A}\widehat{\mathit{x}}}_{k+d+i|k}+{\mathit{B}\widehat{\mathit{u}}}_{k+i}+{\mathit{A}}^{d+i}{\mathit{B}}_e{\mathit{w}}_k,{\widehat{\mathit{y}}}_{k+d+i|k}={\mathit{C}\widehat{\mathit{x}}}_{k+d+i|k},{\widehat{\mathit{x}}}_{d|0}={\mathit{x}}_d,\hfill \\ & & {\widehat{\mathit{y}}}^{min}\le {\widehat{\mathit{y}}}_{k+d+i}\le {\widehat{\mathit{y}}}^{max},{\widehat{\mathit{u}}}^{min}\le {\widehat{\mathit{u}}}_{k+i}\le {\widehat{\mathit{u}}}^{max},\mathsf{\Delta }{\widehat{\mathit{u}}}_{k+i}=0fori\in \{H_u+1,\cdots ,H_p\}.\hfill \end{array}$$

3.6. ESN with Gaussian Process Regression (GPR)

ESN, initially developed by Jaeger, is a machine learning technique based on a recursive neural network [21,22]. This architecture is structured to create a reservoir with a randomly connected recurrent network that encodes the underlying dynamics in the internal states. This configuration enables fast learning as a simple algorithm is used. Furthermore, the training procedure, where only the readout unit is trained, is relatively inexpensive. In the proposed method, ESN was applied to predict future disturbance. The time required for the optimization calculation was ensured to be within the time margin between the current and actual execution times (Figure 2). ESN involves three components, the input units, reservoir, and output units, as shown in Figure 3. The reservoir, which is entirely connected with K input and L output units, transforms the input data into nonlinear features in a high-dimensional space. To obtain the same output as the target value, the output weight matrix was updated in the output layer only. The activations of the input units, internal units, and output units at time step t are $\mathit{s}\left(t\right)$=${(s_1\left(t\right),\cdots ,s_K\left(t\right))}^{\mathrm{T}}$, $\mathit{r}\left(t\right)$=${(r_1\left(t\right),\cdots ,r_N\left(t\right))}^{\mathrm{T}}$, and $\mathit{z}\left(t\right)$=${(z_1\left(t\right),\cdots ,z_L\left(t\right))}^{\mathrm{T}}$, respectively. The internal activations were updated in the reservoir with leaky integrator neurons, as follows:

$$\begin{array}{cc}\hfill \mathit{r}\left(t\right)& =(1-a)\mathit{r}(t-1)+af({\mathit{W}}^{in}\mathit{s}\left(t\right)+\mathit{W}\mathit{r}(t-1))\end{array}$$

(33)

where f is the activation function of the reservoir, and a is the leaking decay rate. The connection weights were randomly generated in an $N\times K$ input weight matrix ${\mathit{W}}^{in}$ and an $N\times N$ internal weight matrix $\mathit{W}$. Generally, these weights are derived using a normal or uniform distribution, which represents fixed values in the network. The output of the reservoir is computed as:

$$\begin{array}{cc}\hfill \mathit{z}\left(t\right)& =f^{out}\left({\mathit{W}}^{out}(\mathit{s}\left(t\right),\mathit{r}\left(t\right),\mathit{z}(t-1))\right)\end{array}$$

(34)

where $f^{out}$ is the activation function of the output neurons. ${\mathit{W}}^{out}$ is an $L\times (K+N+L)$ output weight matrix, which is trained using a linear regression algorithm.

Furthermore, we assumed that the time-series data for learning the ESN included system noise. Additionally, a Gaussian process regression (GPR) [23] was introduced to estimate the model using the input data adaptively. GPR is a method that constructs a regression model of nonlinear time-series data by assuming the joint distribution of each value of time-series by following the Gaussian distribution. Generally, the regression process of time-series data involves estimating the regression model $yi=f\left(t_k\right)$ for arbitrary input $t_k$ when the outputs ${(y_1,y_2,\dots ,y_M)}^{\mathrm{T}}$ for time-series inputs ${(t_1,t_2,\dots ,t_M)}^{\mathrm{T}}$ are given. Assuming that the joint distribution of the M-dimensional vector $\mathit{y}$ consisting of time-series data $\mathit{y}={(y_1,y_2,\dots ,y_m)}^{\mathrm{T}}$ follows a Gaussian distribution. The Gaussian process that provides the models for joint distribution can be described by the mean $E\left[\mathit{y}\right]$ and covariance $Cov\left[\mathit{y}\right]$. By estimating the mean and covariance, a probabilistic model can be generated for the given time-series data. By assuming the observed noise for $\mathit{y}$, the target regression model can be calculated by using the covariance matrix ${\mathit{C}}_M$, whose elements are the sum of the kernel functions $k_f(t_i,t_j)$ that compose the Gram matrix and the observed noise ${\beta }^{-1}{\delta }_{i,j}$ of $\mathit{y}$, where ${\delta }_{i,j}$ is a delta function that becomes 1 when $i=j$. To estimate the output ${\widehat{y}}_{M+1}$ for the new input $t_{M+1}$, we would have to consider the average of the Gaussian distribution at the input $t_{M+1}$. Based on the vector ${\mathit{k}}_f$ and the given time-series data $\mathit{y}$, the predicted mean value ${\widehat{y}}_{M+1}$ and covariance ${\widehat{\sigma }}_{M+1}^2$ were calculated as follows:

$$\begin{array}{ccc}& & {\widehat{y}}_{M+1}={\mathit{k}}_f^{\mathrm{T}}{\mathit{C}}_M^{-1}\mathit{y},\hfill \end{array}$$

(35)

$$\begin{array}{ccc}& & {\widehat{\sigma }}_{M+1}^2=k_f(t_{M+1},t_{M+1})+{\beta }^{-1}-{\mathit{k}}_f^{\mathrm{T}}{\mathit{C}}_M^{-1}{\mathit{k}}_f\hfill \end{array}$$

(36)

where ${\mathit{k}}_f$ has the elements whose kernel function $k_f(t_i,t_{M+1})$ was calculated for each ($i=1,\dots ,M$). We selected the following kernel function:

$$\begin{array}{c}\hfill k_f(t_i,t_j)=a_{h}exp(-b_h{(t_i-t_j)}^2).\end{array}$$

$a_h$ and $b_h$ are the hyperparameters that maximize the logarithmic likelihood, which can be estimated by the gradient descent method. By applying ${\theta }_h={(a_h,b_h)}^{\mathrm{T}}$, the logarithmic likelihood can be written as [23]:

$$\begin{array}{c}\hfill lnp\left(\mathit{y}\right|{\theta }_h)=-\frac{1}{2}ln\left|{\mathit{C}}_M\right|-\frac{1}{2}{\mathit{y}}^{\mathrm{T}}{\mathit{C}}_M^{-1}\mathit{y}-\frac{M}{2}ln2\pi .\end{array}$$

(37)

4. Simulation

To evaluate the control performance, we considered two simulation scenarios involving the comparative models of the proposed ESN-MPC and EDT-MPC. The proposed ESN-MPC model employs the optimization problem formulated in Equation (21). This formulation does not include the effect of the dead-time. The prediction obtained using the ESN provides future disturbance information. To consider the effect of the dead-time, the input signal is applied to the future control. The optimal control for the future is defined at the current time. Instead of solving the optimization problem in real-time, the future optimized signal is buffered and applied at the appropriate time to address the dead-time. The time difference between the current and MPC optimization times, termed the time margin, can be set arbitrarily regardless of the dead-time, as illustrated in Figure 2. The reasonable requirement is that the time margin can be set longer than the sum of the calculation time and dead-time without sacrificing the control performance. In contrast, the conventional EDT-MPC model includes the effect of the dead-time, as indicated in Equation (32). However, the conventional EDT-MPC model does not include the time margin. Instead, to perform a comparison with the proposed method, the performance of the future predictive control was evaluated by setting the artificial dead-time in the EDT-MPC to the same value as the time margin in ESN-MPC.

The performances of the two methods were evaluated via computer simulations. All simulations were performed in MATLAB, R2019a installed on a PC with a clock speed of 3.1 GHz and 16 GB of RAM in a Mac OS 10.12 environment.

4.1. Numerical Conditions

We evaluated the proposed system considering the torque oscillation at the engine start pattern. The waveforms were created using the engine simulator that was provided as a sample model [17]. In real-life operation, the conditions exhibit many variations in terms of engine temperature, humidity, etc. To consider such variations, we set two parameters, i.e., friction torque of the engine and starter speed, to have a normal distribution with a variation of 5 %. The training and test data sets for the ESN were established based on the engine simulation with random noise.

Figure 4a illustrates a sample of typical torque oscillation during the engine start condition. The engine does not burn fuel until approximately 0.8 s. In this no combustion condition, torque oscillation is mainly caused by pumping and gas exchange loss in the engine cylinder. After approximately 0.8 s, the engine torque increases sharply with the initiation of fuel combustion and later decreases to the torque demand for the idle speed. Figure 4b shows the overlap of 50 variations of the waveform. It can be noted that the frequency of the waveform exhibits considerable variations. To evaluate the control performance in cases involving irregular disturbances, we applied the aforementioned waveforms to the active vibration reduction controller. The difference in the control performance between the conventional MPC and the proposed method is discussed in the subsequent section. The specifications of the active vibration controller for the HEV powertrain are presented in

Table 1. Specifications of the active vibration controller for the reference HEV powertrain. The control performances of the proposed ESN-MPC were compared using different design variables. ( ) represents the case where the inertia is 0.08 kg m${}^2$.

Symbol	Parameter	Value
$I_{m1}$	M1 inertia	0.0265 kg m${}^2$
$I_{m2}$	M2 inertia	0.035 kg m${}^2$
$I_e$	Engine inertia	0.20 (0.08) kg m${}^2$
$I_c$	Carrier inertia	0.01 kg m${}^2$
$I_s$	Sun gear inertia	0.01 kg m${}^2$
$I_r$	Ring gear inertia	0.005 kg m${}^2$
$I_p$	Pinion gear inertia	0.001 kg m${}^2$
$R_s$	Sun radius	0.0477 m
$R_p$	Pinion radius	0.0382 m
$n_p$	Number of pinions	4
$\lambda$	Planetary gear ratio	0.3846
$k_{dp}$	Damper stiffness	700 Nm/rad
$c_{dp}$	Damping coefficient	10 Nm s/rad
$t_s$	Sampling time	1 ms
$H_u$	Control horizon	5 steps
$H_p$	Prediction horizon	6 steps
$T_{m1}^{min}$, $T_{m1}^{max}$	M1 torque limitation	−150, +150 Nm
$T_{m2}^{min}$, $T_{m2}^{max}$	M2 torque limitation	−100, +100 (−200, +200) Nm

.

4.2. ESN Prediction Results

The training dataset for the ESN prediction model included 40 realization samples of the engine torque waveform. After learning the dataset, 10 untrained datasets were used for the validation. To verify the application of the ESN, simulation parameters were applied to the model as follows: input unit K$=$1, output unit L$=$1, reservoir units N$=$500, and leaking decay rate a$=$$0.3$. These parameters were determined empirically as the optimal values for the studied benchmarks. The input data includes the time-series waveform. One input unit received the past waveform, and one output unit determined the value for the next time step. By repeating this operation, the trajectory tracking of the waveform could be predicted.

Figure 5 shows the prediction results of the ESN for the trajectory tracking of the waveform. Each figure displays the estimation performed by the ESN for different forward steps from the current time. As expected, the ESN estimation is nearly similar to the actual waveform shown in Figure 5a (forward step in 1 ms). However, as the forward time steps increased to 10 and 20 ms, as shown in Figure 5b,c, the difference between the estimated and true state waveforms increased. To examine this aspect, the total error in one waveform was estimated for each condition. The arithmetic mean of the total error over all points of one waveform was calculated as:

$$\begin{array}{cc}& e=\frac{1}{N_{sample}}{\sum }_{i=1}^{N_{sample}}\frac{\mid T_i^{estimated}-T_i^{actual}\mid }{\mid T_i^{actual}\mid }\end{array}$$

(38)

where $N_{sample}$ is the number of sampling data points in one target waveform, $T_i^{estimated}$ is the estimated result from the ESN, and $T_i^{actual}$ is the target waveform. $T_i^{actual}$ is the output torque from the ICE and includes pumping loss. Pumping loss is energy loss due to the reciprocating motion of the piston. Since part of the energy generated by combustion is consumed in the intake and exhaust stroke, $T_i^{actual}$ takes either a positive or negative value. The results are shown in Figure 6. The error bars show the maximum and minimum data for all 10 runs. The error increased monotonically, thereby indicating that the predictions for the distant future tended to be less accurate.

4.3. ESN-Based MPC Results

After training the ESN, we obtained the trained weight matrix ${\mathit{W}}^{out}$ and applied it to the proposed ESN-MPC scheme. Figure 7 shows the results of the proposed active vibration control for the ICE start condition compared to that for the HEV drivetrain without a controller. The untrained samples were applied to the proposed ESN-MPC scheme. For comparison, the time margin between the current and MPC optimization times (as illustrated in Figure 2) was set as one control period for the proposed controller. Figure 7a,c,e show the results for the proposed controller, whereas Figure 7b,d,f show the results for the case without a controller. Figure 7c,d show the results for the output torque. The proposed controller reduces the torque oscillation compared to that in the case without the controller, despite the occurrence of sudden torque fluctuations with a large amplitude. The rotation speeds ${\omega }_e$ − ${\omega }_c$ and ${\omega }_{m2}$ were nearly stable at 0 rpm during the entire engine starting period, as shown in Figure 7a. Generally, ${\omega }_{m2}$ is related to the output speed, and 0 rpm indicates that the vehicle is not moving. In contrast, in the case without a controller, as shown in Figure 7b, ${\omega }_{m2}$ is stable at 0 rpm, whereas ${\omega }_e$ − ${\omega }_c$ exhibits the typical characteristics of an oscillating waveform. Generally, ${\omega }_e$ − ${\omega }_c$ indicates the difference in the engine rotation and planetary carrier rotation, which corresponds to the angular speed of the torsional damper. Considering the function of the torsional damper, when the input speed is equal to the output speed (${\omega }_e$ − ${\omega }_c$$=$0), the torsional damper does not have a twist angle. In this case, the HEV powertrain does not transmit any torque to drive the output shaft. Therefore, the proposed ESN-MPC controller can reduce the torque oscillation. For the control input, the proposed controller inputs the torque of the two motors, i.e., $T_{m1}$ and $T_{m2}$, as shown in Figure 7e. By comparing Figure 7d,e, it can be noted that the proposed controller tends to generate an inverse torque for the oscillation.

We evaluated EDT-MPC using the conventional approach and compared it with the proposed ESN-MPC. Both the controllers considered the effect of the dead-time and time lag behavior of the actuator in different ways. The time margin applied to the proposed ESN-MPC and the artificial dead-time applied to the conventional EDT-MPC indicates the time difference between the current and actual control times. Generally, this time difference affects the prediction accuracy of the disturbance and compensation performance. Figure 8 shows the results for the proposed ESN-MPC and conventional EDT-MPC during the ICE start condition. Three different time margins (1, 10, and 20 ms) are compared in the results. Figure 8a–c present the output torque of the HEV powertrain. In the 1 and 10 ms cases as shown in Figure 8a,b, the output torque is stable at 0 Nm when the proposed ESN-MPC is used; however, in the case of the conventional EDT-MPC, the torque tends to converge with minor oscillations. In the 20 ms case Figure 8c, the result for the proposed controller is stable, whereas that of EDT-MPC exhibits increased fluctuation, thereby indicating the higher performance of the ESN-MPC controller compared to that of EDT-MPC. Figure 8d–i presents two control inputs. These control inputs remain within the corresponding limits in all the cases owing to the MPC constraints. The control inputs for the proposed method do not exhibit a large difference as the time margin increases, as shown in Figure 8d–f. In contrast, the control inputs when using the conventional EDT-MPC, as shown in Figure 8g–i, exhibit a considerable change as the setting dead-time increases. In particular, as the setting dead-time increases, the amplitude of the control input increases, and the shape of the control input changes. In the 20 ms case, the waveforms of the two control inputs are inverted compared to those at 10 ms, demonstrating the occurrence of a phase shift. The fluctuations of the output torque are larger in the conventional EDT-MPC than those in the proposed method. Owing to the phase shift of the control input and larger amplitude, a larger compensation error occurs while using the conventional controller. The proposed ESN-MPC outperformed the conventional EDT-MPC even when the time margins were large.

The execution time for the optimization computation for the two controllers is compared in

Table 2. The computation time for the optimization of 10 test samples.

			CPU Time (ms)
Model	Simulation Condition		Ave.	Max.
ESN-MPC	Time margin	1 ms	1.2	6.1
		10 ms	1.2	7.7
		20 ms	1.2	5.1
EDT-MPC	Setting dead-time	1 ms	2.2	5.9
		10 ms	2.3	7.5
		20 ms	3.4	9.5

. The histograms of the computation time are shown in Figure 9. The measured computation time includes all the 10 patterns of test data. The proposed ESN-MPC and conventional EDT-MPC were tested for three values of the time margin and setting dead-time, respectively, 1, 10, and 20 ms. In Table 2, the average computation time when using the proposed ESN-MPC is the same, regardless of the time margin. Similarly, the histograms of the computation time shown in Figure 9 exhibit nearly the same shape when using the proposed method. In contrast, the average computation time when using the conventional EDT-MPC increases with the setting dead-time. The distribution of the calculation time when using the conventional method tends to shift towards a larger computation time as the setting dead-time increases. Owing to the dead-time required for the recursive calculation in the conventional EDT-MPC, a larger dead-time corresponds to a larger computation time. Moreover, the maximum computation time is different under different simulation conditions. Although the frequency is extremely low, it is assumed that rare conditions exist in the optimization calculation in which a significant amount of time is required to identify the optimal solution in the test case for both models.

4.4. Case Study for the Design Variables

The predictive model in the proposed method for the HEV system is based on a linear model with the state-space representation to clarify the relationship between the model parameters and outputs. To evaluate the performance and adequacy of the design variables in the proposed method, we compared the results for three cases of the design variables, as shown in Figure 10. All the simulations were performed considering the same time margin of 10 ms in the proposed ESN-MPC. Figure 10a,d,g show the results for the original design as references. The engine inertia was changed to 0.08 kg m${}^2$. The torque limitations of the motor remained constant in (b), (e), and (h). However, the limitations changed in (c), (f), and (i). Decreasing the inertia results in deteriorated control performance, as indicated by the increased torque oscillations observed in Figure 10b. Moreover, the difference in the rotational speed between the engine crankshaft and carrier, ${\omega }_e$ − ${\omega }_c$, begins to oscillate, as shown in Figure 10e, and the motor torque pertaining to the two control inputs attain the limit values, as shown in Figure 10h. In case (c), the torque limits of the motors were higher than those in case (b). This corresponds to an improved vibration control performance. The increased motor limits indicate that the control inputs can implement the appropriate control for vibration reduction.

5. Discussion

In this study, we assumed a realistic system disturbance and verified it using a simulation with noise. The results in Figure 8 indicate that the proposed control scheme based on the combined ESN and MPC can improve the control performance even if future disturbances change drastically under a large time margin. These findings indicate that by using the prediction obtained by the ESN, the proposed method can be applied to predictive control for future timing. However, EDT has been proposed as a method of compensating for the realistic dead-time at the current control timing. In this study, EDT-MPC was extended beyond the original compensation timing for comparison with the proposed method. Hence, the accuracy of compensation for the future timing in EDT-MPC is lower than that of the proposed method. It should be mentioned that the conditions were slightly unfair for comparison. Note that the main focus of this study is to demonstrate that the proposed method can be applied to predictive control for future timing. Moreover, the results of the computation time presented in Figure 9 demonstrate the feasibility of the proposed method. It is a reasonable requirement for the optimized calculation time and dead-time to be smaller than the time margin. The proposed method can be applied when the control performance can be satisfied with an appropriate time margin. As long as the disturbance can be reasonably measured and predicted, this scheme can be applied even if the optimization calculation is not completed within the control period. Moreover, the predictive model in the proposed method for the HEV system employs a linear model with state-space representation. The results shown in Figure 10 suggest that changing the engine inertia changes the torque balance, and appropriate torque limits can be established using the proposed method. In such examinations, physical-based linear modeling can help in conducting a design variable study for numerous cases. However, all models cannot be represented by state-space representations in real systems. Therefore, it is useful to utilize inference-based models using neural networks or similar ESNs.

6. Conclusions

In this study, an ESN-based MPC controller, which calculates the future control input in advance, is proposed to achieve active vibration control of an HEV powertrain. The improvements achieved by the proposed method were discussed under an ICE starting scenario by comparing the findings with those obtained using the conventional EDT compensation scheme by assuming drastic changes in the time-varying disturbances. Different values of the time margin were considered in the simulations. The conventional MPC has a limitation for real-time control, which is that the optimization calculation must be completed within the next control period. In this study, by creating the time margin required for the optimization calculation, we showed the possibility that real-time control can be applied to controlled objects for which optimization cannot be completed within the control period. The simulation results indicated that the proposed method improved the control performance compared to that of the conventional method even when the time margins were large, and the validity of the proposed method was confirmed. Moreover, the proposed method can be applied not only for the considered application but also for other control targets. To examine the system behavior pertaining to the MPC, a linear prediction model based on the state-space representation was used instead of an inference-based model. Disturbance prediction was performed separately using an ESN predictor. Because the model parameters were expressed in the linear system, the relationship between the parameters and their results could be verified. Future research will be focused on enhancing the learning method of time-varying disturbances by considering generalized driving conditions and applications. Furthermore, autonomous vehicles can be expected to consume substantially less gas and energy by eliminating undesirable ways of human driving [24]. Another part of our future work will involve improving comfortable autonomous driving by learning various driving patterns, which is possible to control the situation more finely than before.