Longitudinal Modelling and Control of In-Wheel-Motor Electric Vehicles as Multi-Agent Systems

: This paper deals with longitudinal motion control of electric vehicles (EVs) driven by in-wheel-motors (IWMs). It shows that the IWM-EV is fundamentally a multi-agent system with physical interaction. Three ways to model the IWM-EV are proposed, and each is applicable to certain control objectives. Firstly, a nonlinear model with hierarchical structure is established, and it can be used for passivity-based motion control. Secondly, a linearized model with rank-1 interconnection matrix is presented for stability analysis. Thirdly, a time-varying state-space model is proposed for optimal control using linear quadratic regulator (LQR). The proposed modellings contribute the new understanding of IWM-EV dynamics from the view point of multi-agent-system theory. By choosing the suitable control theory for each model, the complexity level of system design is maintained constant, no matter what the number of IWMs installed to the vehicle body. The e ﬀ ectiveness of three models and their design approaches are discussed by several examples with Matlab / Carsim co-simulator.


Introduction
In-wheel-motors (IWMs) have been shown a novel actuator for motion control system of electric vehicles (EVs) [1]. A remarkable merit of IWM-EV is the ability to precisely generate the driving/braking torques with fast response at individual wheels. Based on this merit, various motion control strategies have been developed in three dimensions. Through literature review, motion control of electric vehicle can be categorized by several sub-fields as follows: (i) Lateral motion control: By differentiating the torques distributed to the left-side and right-side wheels, a yaw-moment is generated about the vertical axis of the vehicle body. This allows us to stabilize the lateral motion of the vehicle via direct-yaw-moment control [2]. Some important issues of this sub-field are sides-slip angle estimation [2,3], side-slip angle and yaw-rate control [4], and the integration of direct-yaw-moment control with active-steering control [5].
(ii) Roll motion control: By allocating the motor torques suitably, it is also possible to generate the roll moment about the longitudinal axis. This moment can be used to control the roll angle and roll rate for improving the comfort of the driver [6].
(iii) Pitch motion control: By distributing the different torques to the front-wheels and rear-wheels, a pitch moment is generated about the lateral axis. This moment could be used to improve the pitch motion, as presented in [7][8][9].  [14][15][16] Optimal control by hierarchical LQR YES [17] (II) Anti-slip control Zero-slip-model following control NO [18,19] Maximum transmissible torque estimation NO [20,21] Direct wheel velocity control NO [22] (III) Driving force control Direct driving force control NO [23] Driving force control based on wheel velocity control and virtual variable control NO [24][25][26] In recent years, an uptrend has been observed in the number of IWM actuators installed to the vehicle body. For instance, France Army developed the armored vehicle DPE 6 × 6 which is driven by six independent motors [27]. In Japan, Keio University introduced the prototype of KAZ-EV with eight IWMs [28]. Thanks to the development of wireless IWM technology [29], powerful EVs driven by more than ten IWMs might be soon realized. With respect to this uptrend, IWM-EV should be treated as a multi-agent system in which each local agent is a locally controlled wheel. The local agent cooperates with each other through the vehicle body to generate the global motion of the vehicle body as well as its local motion.
From the above discussion, LMC of IWM-EVs becomes an actual challenge. Fortunately, it is possible to apply the state-of-the-art of multi-agent system theories. Among them, the global (global/local) framework proposed by Hara group has been shown an effective way to model, analyze, and design multi-agent dynamical systems [30]. Several successful examples of glocal control are gene-regulatory network [31], platoon car [32,33], and manufacturing systems [34]. We recently developed a scheme in the glocal theory family, namely, hierarchical linear quadratic regulator (H-LQR), and applied it to optimal control of slip ratio for IWM-EVs in the acceleration mode [17].
The key of H-LQR is to utilize some special features in the structure of the multi-agent system to reduce the design cost. From the view point of glocal control framework, two critical issues still remain in the research field of LMC for IWM-EVs.
Firstly, almost the previous studies (including the works mentioned in Table 1) neglected the physical interaction between the wheel dynamics. For instance, in [22], Suzuki et al. introduced an anti-slip control by directly controlling the wheel velocity. The wheel velocity controller C w (s) is merely designed by stable pole-placement to the closed-loop of C w (s) and P w (s) which represent the local wheel dynamics. In [25], Maeda et al. proposed the driving force control system for the four-IWM-EV. However, the system is just the combination of four decoupled wheel-velocity-based force control loops. In [11], Nam et al. independently design for each wheel a sliding-mode slip ratio controller. However, the design is based on the single wheel dynamics model without any relevant to the others. Actually, the physical interaction exists as the wheels are installed to the same car body. As clarified by a glocal control tool, namely "generalized frequency variable theory" [35], and its application to EV [36], neglecting physical interaction might reduce the system performance. Even if each wheel's control loop is locally stable, this does not automatically guarantee the stability of the IWM-EV system as a whole. Therefore, it is essential to model the IWM-EV as an interconnected system of multi-agents to properly design the LMC.
Secondly, following the uptrend in the number of IWM actuators, analyzing and designing the LMC systems become more and more complex. Increasing the actuator numbers means increasing the Energies 2020, 13, 5437 3 of 28 order of the system, and hence the implementation effort and the computational cost. Even if we can model the LCM system with physical interactions, how should we select the suitable control theory to reduce the burden of design and analysis?
To tackle the aforementioned issues, proper modelling is necessary. In other words, the model should capture the essential characteristics of the plant in the physical world. Looking at the literature of EV-related-study, this philosophy of modelling has been pointed out by several recent works on battery management and control [37][38][39][40]. For instance, in [37], a predictive control algorithm of charging process was proposed by utilizing a thermoelectric model which accurately captures the battery dynamics. Recently, a Gaussian process regression models have been used to address the aging phenomenon of batteries in [39].
To provide the proper models for LMC of IWM-EVs, this paper is based on the glocal framework. The strategy of this paper is described in Figure 1 which includes three spaces. In the physical world, the motion of the IWM-EV is governed by the physical equations. From the view point of multi-agent systems, the physical model is mapped to the control world by three design models (DMs). Each DM is suitable for a certain objective of longitudinal motion control. Moreover, each DMs is associated with a control theory developed in the glocal framework, such that the complexity of system design and analysis is independent of the number of IWM actuators. Three proposed DMs are briefly explained as follows.
Energies 2020, 13, x FOR PEER REVIEW 3 of 27 Secondly, following the uptrend in the number of IWM actuators, analyzing and designing the LMC systems become more and more complex. Increasing the actuator numbers means increasing the order of the system, and hence the implementation effort and the computational cost. Even if we can model the LCM system with physical interactions, how should we select the suitable control theory to reduce the burden of design and analysis?
To tackle the aforementioned issues, proper modelling is necessary. In other words, the model should capture the essential characteristics of the plant in the physical world. Looking at the literature of EV-related-study, this philosophy of modelling has been pointed out by several recent works on battery management and control [37][38][39][40]. For instance, in [37], a predictive control algorithm of charging process was proposed by utilizing a thermoelectric model which accurately captures the battery dynamics. Recently, a Gaussian process regression models have been used to address the aging phenomenon of batteries in [39].
To provide the proper models for LMC of IWM-EVs, this paper is based on the glocal framework. The strategy of this paper is described in Figure 1 which includes three spaces. In the physical world, the motion of the IWM-EV is governed by the physical equations. From the view point of multi-agent systems, the physical model is mapped to the control world by three design models (DMs). Each DM is suitable for a certain objective of longitudinal motion control. Moreover, each DMs is associated with a control theory developed in the glocal framework, such that the complexity of system design and analysis is independent of the number of IWM actuators. Three proposed DMs are briefly explained as follows.  Figure 1. Discussion strategy in this paper: Design procedure of electric vehicles driven by in-wheelmotors (IWM-EV) motion control.
DM-1 (nonlinear model with hierarchical structure): DM-1 consists of the vehicle body dynamics V in the upper-layer, and a bunch of local wheel dynamics {Wi} in the lower-layer. The two layers interconnect with each other through the aggregation and distribution channels. This expression is quite useful for designing LMC by passivity theory [41]. The advantage of using DM-1 and the passivity theory is that we can design the control system with less mathematical effort and prove the system stability rigorously. In this paper, a passivity based anti-slip control system is proposed to demonstrate the effectiveness of DM-1.
DM-2 (linearized model with rank-1 interconnection matrix): DM-2 is established by linearizing the nonlinear physical equations of IWM-EV about a given operating point, and it is useful for stability analysis. As an example, DM-2 is applied to wheel velocity control with work-load minimization. Thanks to the rank-1 interconnection matrix and the generalized frequency variable (GFV) theory [35], the complexity of stability analysis is considerably reduced. The stability analysis can be performed by a limited set of inequalities established from the GFV transfer function and the eigenvalues of the interconnection matrix. DM-2 and the GFV makes up the useful tools to design and analyze any existing LMC methods.
DM-3 (time-varying state-space model): DM-3 is proposed for capturing the slip ratio dynamics and driving force dynamics of IWM-EV in the state space expression. Hence, this model can be used for applying the optimal control strategies and updating the optimal control gains in real-time operation. In our recent study [17], DM-3 was established in the acceleration mode. To complete the state space modelling for LMC, this paper investigates the DM-3 in the deceleration mode. We theoretically show that global optimality can be achieved by a hierarchical LQR (H-LQR) algorithm DM-1 (nonlinear model with hierarchical structure): DM-1 consists of the vehicle body dynamics V in the upper-layer, and a bunch of local wheel dynamics {W i } in the lower-layer. The two layers interconnect with each other through the aggregation and distribution channels. This expression is quite useful for designing LMC by passivity theory [41]. The advantage of using DM-1 and the passivity theory is that we can design the control system with less mathematical effort and prove the system stability rigorously. In this paper, a passivity based anti-slip control system is proposed to demonstrate the effectiveness of DM-1.
DM-2 (linearized model with rank-1 interconnection matrix): DM-2 is established by linearizing the nonlinear physical equations of IWM-EV about a given operating point, and it is useful for stability analysis. As an example, DM-2 is applied to wheel velocity control with work-load minimization. Thanks to the rank-1 interconnection matrix and the generalized frequency variable (GFV) theory [35], the complexity of stability analysis is considerably reduced. The stability analysis can be performed by a limited set of inequalities established from the GFV transfer function and the eigenvalues of the interconnection matrix. DM-2 and the GFV makes up the useful tools to design and analyze any existing LMC methods.
DM-3 (time-varying state-space model): DM-3 is proposed for capturing the slip ratio dynamics and driving force dynamics of IWM-EV in the state space expression. Hence, this model can be used for applying the optimal control strategies and updating the optimal control gains in real-time operation. In our recent study [17], DM-3 was established in the acceleration mode. To complete the state space modelling for LMC, this paper investigates the DM-3 in the deceleration mode. We theoretically show that global optimality can be achieved by a hierarchical LQR (H-LQR) algorithm in which we only Energies 2020, 13, 5437 4 of 28 need to solve the local Riccati equation. Applying the proposed H-LQR, the computational cost is effectively reduced.
The remainder of this paper is organized as follows. Section 2 presents the physical equations that describing the longitudinal motion of IWM-EV. To clarify the motivation of this paper, Section 3 discusses the disadvantages of the traditional LMC methods by numerical simulations. Then, DM-1, DM-2, and DM-3 are investigated in Sections 4-6, respectively. Each of the proposed DMs is established under reasonable assumptions, following by a motion control example. Section 7 presents the simulation results of the motion control examples in Sections 4-6. The simulations were conducted using a Matlab/Casim co-simulator which is based on an experimental electric vehicle developed by our research group. Finally, the concluding remarks are given in Section 8.

Longitudinal Dynamics of Vehicle
This paper examines the longitudinal motion of IWM-EV illustrated in Figure 2 of which the number N of IWMs is an even number. Let m be the vehicle mass. Under the assumption that the wheels are homogeneous, J w and r are the inertia and radius of the wheel, respectively. T i , F i , and Z i are the motor torque, the driving force, and the vertical force acting at the ith wheel, respectively. The local motion of the ith wheel is represented by the rotational velocity ω i . The global motion of the vehicle body is represented by the longitudinal velocity v x . F air is the air resistance acting on the vehicle body. The total influence of gravity and rolling resistance is represented by F ξ . that describing the longitudinal motion of IWM-EV. To clarify the motivation of this paper, Section 3 discusses the disadvantages of the traditional LMC methods by numerical simulations. Then, DM-1, DM-2, and DM-3 are investigated in Sections 4-6, respectively. Each of the proposed DMs is established under reasonable assumptions, following by a motion control example. Section 7 presents the simulation results of the motion control examples in Sections 4-6. The simulations were conducted using a Matlab/Casim co-simulator which is based on an experimental electric vehicle developed by our research group. Finally, the concluding remarks are given in Section 8.

Longitudinal Dynamics of Vehicle
This paper examines the longitudinal motion of IWM-EV illustrated in Figure 2 of which the number N of IWMs is an even number. Let m be the vehicle mass. Under the assumption that the wheels are homogeneous, Jw and r are the inertia and radius of the wheel, respectively. Ti, Fi, and Zi are the motor torque, the driving force, and the vertical force acting at the ith wheel, respectively. The local motion of the ith wheel is represented by the rotational velocity ωi. The global motion of the vehicle body is represented by the longitudinal velocity vx. Fair is the air resistance acting on the vehicle body. The total influence of gravity and rolling resistance is represented by Fξ.
The motion of the vehicle body is given by Neglecting the wind velocity, the air resistance is expressed as where the coefficient ρ is calculated from the density of air, the cross-sectional area of the vehicle in the air, and the drag coefficient which depends on the shape of the vehicle [42]. The rotational motion of the ith wheel is given as Let ε be a small number to avoid division-by-zero, the slip ratio of the ith wheel is defined as The nonlinear relationship between the driving force and the slip ratio can be described by several empirical functions. The most popular one is the magic formula proposed by Pacejka [43]. This formula is given by The motion of the vehicle body is given by Neglecting the wind velocity, the air resistance is expressed as where the coefficient ρ is calculated from the density of air, the cross-sectional area of the vehicle in the air, and the drag coefficient which depends on the shape of the vehicle [42]. The rotational motion of the ith wheel is given as Let ε be a small number to avoid division-by-zero, the slip ratio of the ith wheel is defined as Energies 2020, 13, 5437

of 28
The nonlinear relationship between the driving force and the slip ratio can be described by several empirical functions. The most popular one is the magic formula proposed by Pacejka [43]. This formula is given by where A i , B i , C i , and D i are the shape factors, with A i = µ i Z i where µ i is the friction coefficient. Note that the friction coefficient and the shape factors are all identifiable from experimental data. Now, the driving force is expressed as The typical curve of the driving force is shown in Figure 3. The driving force gradually increases as the slip ratio increases from 0 to an optimal value λ* which depends on the road condition. Since the slip ratio exceeds λ*, the driving force is decreased to a minimum value as the slip ratio approach to 1. Investigating the magic formula, we suggest the following practical remark.
where Ai, Bi, Ci, and Di are the shape factors, with Ai = μiZi where μi is the friction coefficient. Note that the friction coefficient and the shape factors are all identifiable from experimental data. Now, the driving force is expressed as The typical curve of the driving force is shown in Figure 3. The driving force gradually increases as the slip ratio increases from 0 to an optimal value λ* which depends on the road condition. Since the slip ratio exceeds λ*, the driving force is decreased to a minimum value as the slip ratio approach to 1. Investigating the magic formula, we suggest the following practical remark.

Motivating Example
In our recent study [17], we examined the slip ratio control method proposed in [14] by both theoretical discussion and numerical simulations. According to [14], it is possible to linearize the dynamics from the motor torque Ti to the slip ratio λi. Then, a PI controller can be designed by placing the stable poles to the local subsystem. However, the performance of the overall system should be determined by both local subsystem's transfer function and a matrix that representing physical interaction [35]. Due to the change of operating points, the poles of the overall system might move towards the imaginary axis with increasing imaginary parts. This phenomenon might degrade the system performance by several scenarios. For instance, oscillated phenomena can be observed when the vehicle runs at high velocity. In addition, the load transfer can introduce large gap between the slip ratios of the front-wheels and rear-wheels. Moreover, the traditional slip ratio control system is highly sensitive to the non-flatness of the road surfaces.
This paper neglects to re-simulate the aforementioned slip ratio control method. Instead, we examine another LMC method, namely, anti-slip control by wheel velocity control [22]. To this end, this paper utilizes Carsim-a standard software for simulating the vehicle dynamics [44]. By using Carsim, it is convenient to perform the control algorithms under different environmental conditions, including the dangerous conditions. The simulation model is built based on the four-wheel experimental electric vehicle developed from the Mitsubishi's i-MiEV prototype. It was an outcome of our research project KC.03.08/11-15 supported by the Ministry of Science and Technology of Vietnam. The main specification of the vehicle is summarized in Appendix A.
To establish the simulation model, we modify the parameters of a hatchback car provided in the Library of Carsim. Since the Carsim only provide the vehicle of gasoline type, we constructed the electrical drive system in Matlab/Simulink. The electrical drive system consists of the battery model, the inverter models, the in-wheel-motor models of interior permanent magnet type, together with the inverter control units and motor control units. To resemble the behavior of the actual drive system,

Motivating Example
In our recent study [17], we examined the slip ratio control method proposed in [14] by both theoretical discussion and numerical simulations. According to [14], it is possible to linearize the dynamics from the motor torque T i to the slip ratio λ i . Then, a PI controller can be designed by placing the stable poles to the local subsystem. However, the performance of the overall system should be determined by both local subsystem's transfer function and a matrix that representing physical interaction [35]. Due to the change of operating points, the poles of the overall system might move towards the imaginary axis with increasing imaginary parts. This phenomenon might degrade the system performance by several scenarios. For instance, oscillated phenomena can be observed when the vehicle runs at high velocity. In addition, the load transfer can introduce large gap between the slip ratios of the front-wheels and rear-wheels. Moreover, the traditional slip ratio control system is highly sensitive to the non-flatness of the road surfaces.
This paper neglects to re-simulate the aforementioned slip ratio control method. Instead, we examine another LMC method, namely, anti-slip control by wheel velocity control [22]. To this end, this paper utilizes Carsim-a standard software for simulating the vehicle dynamics [44]. By using Carsim, it is convenient to perform the control algorithms under different environmental conditions, including the dangerous conditions. The simulation model is built based on the four-wheel experimental electric vehicle developed from the Mitsubishi's i-MiEV prototype. It was an outcome of our research project KC.03.08/11-15 supported by the Ministry of Science and Technology of Vietnam. The main specification of the vehicle is summarized in Appendix A.
Energies 2020, 13, 5437 6 of 28 To establish the simulation model, we modify the parameters of a hatchback car provided in the Library of Carsim. Since the Carsim only provide the vehicle of gasoline type, we constructed the electrical drive system in Matlab/Simulink. The electrical drive system consists of the battery model, the inverter models, the in-wheel-motor models of interior permanent magnet type, together with the inverter control units and motor control units. To resemble the behavior of the actual drive system, the inverter model imitates the electrical circuit of the experimental inverter. The output torques of the drive system model in Matlab/Simulink are sent to Carsim through the imports namely IMP_MYUSM_*. The notation * represents the IWM position (L1, L2, R1, and R2, where L denote "left" and R denotes "right").

Anti-Slip Control Based on Wheel Velocity Control
Following the idea of the original work [22], we generalize the block diagram of the control system as in Figure 4. Here, 1 N is the all-one column-vector of size N, and The scalar signal T cmd is the driving command given by the driver or the upper motion control layer. In addition, vector k = [k 1 . . . k N ] T is to distribute the driving command to the local wheel. The distribution ratio k i are commonly selected such that 0 < k i < 1 and k 1 + k 2 + . . . + k N = 1. We assume that the vehicle operates on the acceleration mode when rω i > v x , and let λ * be the desired slip ratio. From (4), the reference velocity can be calculated from the wheel velocity v x via the gain the inverter model imitates the electrical circuit of the experimental inverter. The output torques of the drive system model in Matlab/Simulink are sent to Carsim through the imports namely IMP_MYUSM_*. The notation * represents the IWM position (L1, L2, R1, and R2, where L denote "left" and R denotes "right").

Anti-Slip Control Based on Wheel Velocity Control [22]
Following the idea of the original work [22], we generalize the block diagram of the control system as in Figure 4. Here, 1N is the all-one column-vector of size N, and T = [T1 … TN] T , ω = [ω1 … ωN] T . The scalar signal Tcmd is the driving command given by the driver or the upper motion control layer. In addition, vector k = [k1 … kN] T is to distribute the driving command to the local wheel. The distribution ratio ki are commonly selected such that 0 < ki < 1 and k1 + k2 + … + kN = 1. We assume that the vehicle operates on the acceleration mode when ω > i x r v , and let λ * be the desired slip ratio.
From (4), the reference velocity can be calculated from the wheel velocity vx via the gain ( ) Each wheel is provided a velocity controller Cw(s). Let Pw(s) = 1/(Jws) be transfer function that characterizes the local dynamics of the wheel. To design Cw(s), the traditional way is to assign the stable poles to the closed loop transfer function For instance, let κ 1 and κ 2 be the stable poles of Pcl(s), Cw(s) becomes a PI controller with the P gain and I gain calculated as ( )

Simulation and Discussion
To evaluate the above design strategy, a test is conducted as follows. A constant driving command Tcmd is equally distributed to four wheels such that the EV accelerates from the initial position on the high friction surface (μ = 0.8). It enters the low friction surface (μ = 0.2) at about 4 s. Two test cases are performed, and their results are summarized in Figure 5. Each wheel is provided a velocity controller C w (s). Let P w (s) = 1/(J w s) be transfer function that characterizes the local dynamics of the wheel. To design C w (s), the traditional way is to assign the stable poles to the closed loop transfer function For instance, let κ 1 and κ 2 be the stable poles of P cl (s), C w (s) becomes a PI controller with the P gain and I gain calculated as Energies 2020, 13, 5437 7 of 28

Simulation and Discussion
To evaluate the above design strategy, a test is conducted as follows. A constant driving command T cmd is equally distributed to four wheels such that the EV accelerates from the initial position on the high friction surface (µ = 0.8). It enters the low friction surface (µ = 0.2) at about 4 s. Two test cases are performed, and their results are summarized in Figure 5. The simulation results in this sub-section suggest that local stabilization does not guarantee global stability of the control system. By using only the local wheel dynamics to design the controller, we cannot explain such performance of the vehicle in Case B. Therefore, it is necessary to model and design the motion control system of IWM-EV as a whole. Three models that effectively capture the physical interaction will be presented in the following three sections.

DM-1: Nonlinear Model with Hierarchical Structure
This section proposes a way to directly model the IWM-EV system using its nonlinear equations. To design the control system, our strategy is based on the passivity theory which has been widely utilized in various fields of control engineering for three decades. As an example, DM-1 is utilized to design an anti-slip control system. Notice that, it is possible to neglect the disturbance Fξ in the theoretical discussion.

Hierarchical Model of IWM-EV
Based on (3)-(6), the model of the local agent Wi is shown as in Figure 6a. It includes two inputs {vx, Ti} and two outputs {Fi, ωi}. The vehicle body dynamics V is obtained from (1) and (2). V has the input  N 1F which is the aggregation of the local driving force, and the output vx which is distributed to the local Wi. Combining {Wi} in the lower-layer and V in the upper-layer, the DM-1 is established as in Figure 6b with the hierarchical structure.
Model of a single wheel dynamics. Case A: Two stable poles of P cl (s) are placed at κ 1,2 = −10 ± 2j (Figure 5a). The wheel velocities are maintained closed to the vehicle velocity. This means the anti-slip function is successfully attained in this case.
Case B: Two stable poles of P cl (s) are placed at κ 1,2 = −100 ± 0j (Figure 5b). The overall system becomes unstable when the vehicle operates on the low-friction surface.
The simulation results in this sub-section suggest that local stabilization does not guarantee global stability of the control system. By using only the local wheel dynamics to design the controller, we cannot explain such performance of the vehicle in Case B. Therefore, it is necessary to model and design the motion control system of IWM-EV as a whole. Three models that effectively capture the physical interaction will be presented in the following three sections.

DM-1: Nonlinear Model with Hierarchical Structure
This section proposes a way to directly model the IWM-EV system using its nonlinear equations. To design the control system, our strategy is based on the passivity theory which has been widely utilized in various fields of control engineering for three decades. As an example, DM-1 is utilized to design an anti-slip control system. Notice that, it is possible to neglect the disturbance F ξ in the theoretical discussion.

Hierarchical Model of IWM-EV
Based on (3)-(6), the model of the local agent W i is shown as in Figure 6a. It includes two inputs {v x , T i } and two outputs {F i , ω i }. The vehicle body dynamics V is obtained from (1) and (2). V has the input 1 T N F which is the aggregation of the local driving force, and the output v x which is distributed to the local W i . Combining {W i } in the lower-layer and V in the upper-layer, the DM-1 is established as in Figure 6b with the hierarchical structure.
Energies 2020, 13, 5437 8 of 28 Based on (3)-(6), the model of the local agent Wi is shown as in Figure 6a. It includes two inputs {vx, Ti} and two outputs {Fi, ωi}. The vehicle body dynamics V is obtained from (1) and (2). V has the input Τ N 1 F which is the aggregation of the local driving force, and the output vx which is distributed to the local Wi. Combining {Wi} in the lower-layer and V in the upper-layer, the DM-1 is established as in Figure 6b with the hierarchical structure.

Passivity Theory and Its Application to DM-1
Consider a system given by a state space equation with the input vector u(t) ∈ ℝ p , output vector y ∈ ℝ p , and the state vector x(t) ∈ ℝ n . Proposition 3: The IWM-EV system in Figure 6b is passive from the input T to the output ω.

Proof of Proposition 3: Define the storage function
then with respect to Remark 1, it is transparent that the inequality S Τ ≤ ω T  holds true. The detailed calculation is neglected for the sake of paper length.
Passivity-based-control (PBC) was introduced three decades ago to the adaptive control problem of rigid robot [45]. It has been widely applied in various fields of control engineering, especially in

Passivity Theory and Its Application to DM-1
Consider a system given by a state space equation with the input vector u(t) ∈ R p , output vector y ∈ R p , and the state vector x(t) ∈ R n .

Definition 1 [41]:
The system (10) is said to be passive from input u to output y if there exists a positive semidefinite function S: R n → R + , called storage function, such that the inequality . S ≤ y T u holds for all x ∈ R n and u ∈ R p . In addition, (10) is input strictly passive (ISP) if . S ≤ y T u − δ u u 2 for some δ u > 0, and output strictly passive (OSP) if . S ≤ y T u − δ y y 2 for some δ y > 0.

Proposition 1:
The IWM-EV system in Figure 6b is passive from the input T to the output ω.

Proof of Proposition 1: Define the storage function
then with respect to Remark 1, it is transparent that the inequality . S ≤ ω T T holds true. The detailed calculation is neglected for the sake of paper length. Passivity-based-control (PBC) was introduced three decades ago to the adaptive control problem of rigid robot [45]. It has been widely applied in various fields of control engineering, especially in multi-agent systems [46]. Following the fundamental passivity theorems [41], if the plant is passive, system stability can be guaranteed by connecting the plant with a strictly passive feedback controller. Motivated by the Proposition 1, there exists many ways to design motion control of IWM-EV using PBC [47]. A design example of anti-slip control will be examined in the next subsection.

Example of DM-1: Passivty Based Anti-Slip Control of IWM-EV
This sub-section demonstrates anti-slip control as an application of the DM-1 and passivity theory. The significance and necessity of anti-slip control can be explained from two points of view as follows.
From a safety point of view: As discussed by Hori et al. in [14], not only the driving force but also the side force strongly depends on the slip ratio. The side force becomes smaller as the slip ratio becomes bigger. If the slip ratio increases considerably due to the change of road condition, the side force becomes drastically smaller. This causes serious problems, including drift-out in front-wheel-drive vehicles, spin in rear-wheel-drive vehicles, and drift-out with rotation in four-wheel-drive vehicles. Such situations are extremely dangerous for the drivers.
From an energy point of view: For simplicity in discussion, we consider a single-wheel vehicle model. Following the Proposition 1, its energy storage function is With respect to the definition of slip ratio expressed in (4), we obtain In the acceleration mode (rω > v x ), a sharp increasing of λ means a sharp increasing of S. Thus, to attain a given velocity pattern v x in the acceleration mode, it is required to use more electric energy from the battery. Next, we examine the deceleration mode (rω < v x ). In this case the side of λ is negative. A drastic decreasing of λ (or increasing of −λ) results in a drastic decreasing of S. Thus, to attain a given velocity pattern v x in the deceleration mode, only small amount of energy can be regenerated back to the battery. In summary, the driving range of electric vehicles would be considerably reduced if we wasted energy by neglecting anti-slip control. The proposed control system is shown in Figure 7 where the scalar T cmd is the driving command given by the driver. The command distributed to each W i is T r,i = k i × T cmd where the set of distribution ratio satisfies k i > 0 for all i from 1 to N, and k 1 + k 2 + ... + k N = 1. Consequently, vectors T r and k are defined as T r = [T r,1 T r,2 ... T r,N ] T and k = [k 1 k 2 ... k N ] T , respectively. Note that vector k can be updated in real-time to attain the additional control objective, such as yaw moment control or minimization of energy consumption. Let K a and K ω be the control gains. Considering the acceleration mode, a possible option of the local anti-slip control law C l,i is given by Energies 2020, 13, 5437 10 of 28 ω ω . Following the Definition 2, the system in Figure 7 is shown to be OSP. This completes the proof.

DM-2: Linearized Model with Rank-1 Interconnection Matrix
The previous section focuses on the design problem using DM-1. Thanks to the passivity theory, we can rigorously prove the stability of the control system even if it is nonlinear and quite complex. However, passivity theory limits the class of controllers that we can select. Therefore, it is also essential to consider the problem of stability analysis. In other words, given a motion control system, how should we evaluate its stability? This section introduces a design model which can reduce the complexity level of the analysis problem. As an example, it is applied to a wheel velocity control system of IWM-EV. Note that, the air resistance and the other disturbances are neglected in the discussion in this section. For the sake of paper length, we only examine the vehicle motion in the acceleration mode, such that max{rωi, vx, ε} = rωi. We focus our stability analysis to the operating point O = {ωo,1, ωo,2, …, ωo,N, vxo} with the following assumption.

Linearized Model of IWM-EV
Under the Assumption 5, the driving force can be expressed as The stability of the control system is guaranteed by the following theoretical result.

Proposition 2:
Considering the acceleration mode, the anti-slip control system in Figure 7 with the local controller C l,i given by (14) is OSP from T r to ω if the control gains are selected such that K a > 0 and K ω > 0.

Proof of Proposition 2:
With respect to the storage function (11), the control law (14), and the dynamics of W i described by (3)-(6), we have There exist four terms on the right-hand side of (15). It is clear from Remark 1 that the first term is non-negative. Similarly, the second term is non-negative. With the positive control gains, Following the Definition 1, the system in Figure 7 is shown to be OSP. This completes the proof.

DM-2: Linearized Model with Rank-1 Interconnection Matrix
The previous section focuses on the design problem using DM-1. Thanks to the passivity theory, we can rigorously prove the stability of the control system even if it is nonlinear and quite complex. However, passivity theory limits the class of controllers that we can select. Therefore, it is also essential to consider the problem of stability analysis. In other words, given a motion control system, how should we evaluate its stability? This section introduces a design model which can reduce the complexity level of the analysis problem. As an example, it is applied to a wheel velocity control system of IWM-EV. Note that, the air resistance and the other disturbances are neglected in the discussion in this section. For the sake of paper length, we only examine the vehicle motion in the acceleration mode, such that max{rω i , v x , ε} = rω i . We focus our stability analysis to the operating point O = {ω o,1 , ω o,2 , . . . , ω o,N , v xo } with the following assumption.
Assumption 1: (i) The vehicle operates on the uniform road condition and the friction coefficients of all wheels are the same. (ii) The slip ratio is small, and the driving force is linearized as F i ≈ S n λ i where S n is the nominal driving stiffness. At O, the values of max{rω o,i , v xo , ε} are almost the same. In other words, we might approximate

Linearized Model of IWM-EV
Under the Assumption 1, the driving force can be expressed as where S n ≈ κ n S n , and η i = rω i − v x . Based on (17), the linearized model of IWM-EV (DM-2) is established as in Figure 8 where η = [η 1 . . . η N ] T and Energies 2020, 13, x FOR PEER REVIEW 11 of 27 Remark 6: The interconnection matrix A is a matrix of rank-1 with N-1 zero eigenvalues, and only one non-zero eigenvalue which equals to -N.

Example of DM-2: Wheel Velocity Control System and Its Stability Condition
This section considers a wheel velocity control system for distributing the driving command Tcmd with a set of distribution ratios {ki}. From (3), the reference of wheel acceleration is calculated as [48] where ˆi F is the obtained by the driving force observer (DFO) as follows where the low-pass-filter is expressed as Q(s) = 1/(τfs + 1) with the small time constant τf. From (19)- (21) and the DM-2, block diagram of the wheel velocity control system is shown in Figure 9 where Cw(s) is the wheel velocity controller. Notice that we neglect to discuss the algorithms to optimize in this section. Such algorithms can be founded in several other works on range extension control of IWM-EV [48][49][50]. The main goal of this section is to provide a suitable model for stability analysis. About the operating point O, the stability of the system in Figure 9 can be discussed from the system in Figure 10 where the transfer function from y to z is diagonalized by

Remark 2:
The interconnection matrix A is a matrix of rank-1 with N-1 zero eigenvalues, and only one non-zero eigenvalue which equals to -N.

Example of DM-2: Wheel Velocity Control System and Its Stability Condition
This section considers a wheel velocity control system for distributing the driving command T cmd with a set of distribution ratios {k i }. From (3), the reference of wheel acceleration is calculated as [48] .
Energies 2020, 13, 5437 12 of 28 whereF i is the obtained by the driving force observer (DFO) as followŝ where the low-pass-filter is expressed as Q(s) = 1/(τ f s + 1) with the small time constant τ f . From (19)- (21) and the DM-2, block diagram of the wheel velocity control system is shown in Figure 9 where C w (s) is the wheel velocity controller. Notice that we neglect to discuss the algorithms to optimize in this section. Such algorithms can be founded in several other works on range extension control of IWM-EV [48][49][50].
The main goal of this section is to provide a suitable model for stability analysis. About the operating point O, the stability of the system in Figure 9 can be discussed from the system in Figure 10 where the transfer function from y to z is diagonalized by

Stability Test
To demonstrate the Stability condition 7, this section considers the PI controller Substitute (24) to (22), the transfer function of the GFV is derived as φ + + + + = = + + +

Stability Test
To demonstrate the Stability condition 7, this section considers the PI controller Substitute (24) to (22), the transfer function of the GFV is derived as φ + + + + = = + + +   Figure 10. Representation of wheel velocity control system using DM-2 as a multi-agent-system.
The system in Figure 10 is a multi-agent-system of N agent H(s), interconnecting via the matrix A.

The system is stable at O if and only if
for all eigenvalues ν i of matrix A. We notice that the zero eigenvalues certainly satisfy (23). Following the generalized frequency variable theory (GFV) [35] with respect to Remark 2, (23) can be checked by the following condition.

Stability condition 1: The wheel velocity control system is stable at the operating point O if the point (-N, 0)
is located in the stable domain Ω c + defined as Ω + := φ(C + ), Ω c + := C\Ω + where the GFV is φ(s) = H −1 (s).

Stability Test
To demonstrate the Stability condition 1, this section considers the PI controller C w (s) = K P s + K I s .
Substitute (24) to (22), the transfer function of the GFV is derived as where the set of {a i , b i } is expressed in the Appendix B. Define the polynomial Following the Theorem 1 in [35], a set of {D k > 0} that defining the stable domain is obtained as where Based on the GFV theory [35], the Stability condition 1 can be tested as follows.

Stability test 1: We consider an operating point O
We establish the linearized model DM-2 under the Assumption 1. We calculate the GFV using (25), and then establish the set of {D k } using (26) and (27). The wheel velocity control system in Figure 9 is stable at O if The complexity level of the Stability test 1 is independent of the number of IWM actuators. At each operating point, the test is performed by a limited set of inequalities (28).

Discussion: Possibility of Robust Stability Test
By using Assumption 1, the LMC system of electric vehicle is finally modelled as a homogeneous multi-agent system as in Figure 10. However, in actual operation, the friction coefficients of all wheel might be not the same. In addition, the friction coefficients might vary due to the change of road condition. Considering this scenario, we can express the LMC system as a multi-agent system with perturbation as in Figure 11a. The perturbation matrix is ∆(s) = diag{∆ i (s)}. We can assume that each local perturbation satisfies ∆ i (s) ∞ ≤ ξ where the volume ξ of perturbation represents the heterogeneous level between the agents and the uncertainties due to the change of road condition. With respect to the structure of matrix A given by (18), there exists the unitary matrix U that diagonalizes matrix A. Thus, we can transform the system in Figure 11a to the equivalent system in Figure 11b where {λ i } is the set of eigenvalues of matrix A. The transfer function from w to z is diagonalized by N transfer functions Energies 2020, 13, x FOR PEER REVIEW 14 of 27 discussion for future study, as the goal of this paper is only to discuss several ways to model the vehicle dynamics such that longitudinal motion control can be designed properly.

DM-3: Time-Varying State-Space Model
In Section 4, the DM-1 is established based on the "physical equations" of the IWM-EV. Section 5 proposes the DM-2 expressed by "transfer functions and interconnection matrix". This section shows that the dynamics of IWM-EV can also be captured by a "state-space model". We examine the vehicle motion in the deceleration mode, such that max{rωi, vx, ε} = vx. Notice that it is possible to neglect the unknown disturbances of Fair and Fξ in the expression of this model.

State-Space Model of IWM-EV in the Deceleration Mode
Calculating the derivative of the slip ratio (4) with respect to (1) and (3), we have Figure 11. Vehicle control system as a multi-agent-system with perturbation.
As ∆ i (s) ∞ ≤ ξ, it follows that ∆(s) ∞ ≤ ξ , and therefore ∆ (s) ∞ ≤ ξ since U is a unitary matrix. Applying the small gain theorem, the perturbed system is robust stable if The zero eigenvalues certainly satisfy (30). With respect to Remark 2, the robust stability condition is reduced to NH n (s) In summary, the DM-2 is still applicable to deal with the uncertainties. In this situation, we have to design for each wheel a local controller such that a triad of conditions are satisfied: (i) the stability condition of the nominal system given by (28); (ii) the model matching condition given by the norm inequality ∆ i (s) ∞ ≤ ξ; and (iii) the robust stability condition given by (31).
Another way to deal with robust stability is to utilize "sector bounded nonlinearity" of the tire-force characteristics ( Figure 12). We can directly apply the Circle and Popov criteria which was presented in Khalil's famous book on nonlinear control [41]. Fortunately, the special structure of matrix A could be utilized to reduce the complexity of "Circle and Popov criteria". We let this discussion for future study, as the goal of this paper is only to discuss several ways to model the vehicle dynamics such that longitudinal motion control can be designed properly.

DM-3: Time-Varying State-Space Model
In Section 4, the DM-1 is established based on the "physical equations" of the IWM-EV. Section 5 proposes the DM-2 expressed by "transfer functions and interconnection matrix". This section shows that the dynamics of IWM-EV can also be captured by a "state-space model". We examine the vehicle motion in the deceleration mode, such that max{rωi, vx, ε} = vx. Notice that it is possible to neglect the unknown disturbances of Fair and Fξ in the expression of this model.

State-Space Model of IWM-EV in the Deceleration Mode
Calculating the derivative of the slip ratio (4) with respect to (1) and (3), we have (32) Figure 12. Sector bounded nonlinearities to capture the uncertainty of tire-force characteristics.

DM-3: Time-Varying State-Space Model
In Section 4, the DM-1 is established based on the "physical equations" of the IWM-EV. Section 5 proposes the DM-2 expressed by "transfer functions and interconnection matrix". This section shows that the dynamics of IWM-EV can also be captured by a "state-space model". We examine the vehicle motion in the deceleration mode, such that max{rω i , v x , ε} = v x . Notice that it is possible to neglect the unknown disturbances of F air and F ξ in the expression of this model.

State-Space Model of IWM-EV in the Deceleration Mode
Calculating the derivative of the slip ratio (4) with respect to (1) and (3), we have For the sake of estimation and control design, it is possible to use the first order dynamics model [51] of the driving force represented by where τ i is the relaxation time constant and S i is the driving stiffness which is a parameter for capturing the linear region of the Pacejka model. Here, S i can be identified from IWM's torque and on-board sensor [52]. Selecting the slip ratio and driving force as state variables, and considering wheel velocities as time-varying parameter, from (30) and (31), the dynamical equation of the ith agent is expressed as follows where Define , the DM-3 is expressed as Energies 2020, 13, 5437 The full block matrix A G (t) represents the physical interaction between the local agents. The acceleration-mode-version of DM-3 when max{rω i , v x , ε} = rω i can be found in our recent study [17]. DM-3 can be used for controlling either driving force or slip ratio. To this end, we proposed in [17] a hierarchically decentralized LQR algorithm which was verified by real-time experiments using an EV driven by two IWMs.

Example of DM-3: Slip Ratio Control in the Deceleration Mode
To demonstrate the practical application of DM-3, we investigate the slip ratio control in the deceleration mode. The following assumptions are considered for the sake of simplicity in presenting our main idea. Assumption 2 is reasonable since the driving force can be estimated accurately using DFO presented in the previous section. From DFO and the estimated/calculated slip ratio, the driving stiffness can be identified off-line or on-line [52]. On the other hand, the vehicle velocity and wheel slip ratio can be estimated using on-board sensor, such as GPS receiver, accelerometer, and wheel encoders. Several methods for estimating slip ratio can be found in [17,53]. We neglect to focus on state estimation and parameter identification, since it is not the main goal of this paper.
The parameters {τ i , S i } depend on the road condition at each wheel and load transfer. This means the parameters {τ i , S i } are not always the same in real-time operation. The Assumption 3 is reasonable if the wheels are symmetrically distributed in the vehicle body, and the local road conditions of the wheels are quite uniform. Since we implement the feedback control using DM-3, the parameter difference among the local agents can be compensated.
For controlling the slip ratios such that they follow the reference value λ * i , the tracking errors are introduced as the augmented states. In the following, time stamping "t" is eliminated for simplifying the expression. The augmented dynamical model of the ith agent is given as State-space modelling of the total system is expressed as follows using Kronecker product ⊗: The LQR controller of the augmented system (39) can be obtained by minimizing the quadratic cost function where t 0 is the initial time, t f is the finished time. Q and S are symmetric positive semidefinite matrices, and R is strictly symmetric positive definite matrix. Q, H, and R are used to penalize the transient state deviation, the final state, and the control effort, respectively. By solving the Riccati equation that satisfies the boundary condition the control gain is obtained as Equation (41) with condition (42) can be solved using the Hamiltonian established from the system model [54]. To reduce the computational effort, especially when N is a big number, we recently proposed in [17] that the feedback gain can be calculated as by appropriate choices of the weighting matrices. As can be seen in Figure 13, the proposed LQR controller has the hierarchical structure. In the lower layer, the local feedback gain K 1 is used to improve the local control performance of slip ratio. The upper-layer has two term, K g1 and K g2 . K g1 is to deal with the physical interconnection represented by Γ N . On the other hand, K g2 is introduced to improve some other global performances, such as the consensus of slip ratios. Since K 1 is obtained by solving the Riccati equation of a local system, the computation cost is independent of the number of IWM actuators. See Appendix C for the details of the hierarchical LQR algorithm.
Energies 2020, 13, x FOR PEER REVIEW 17 of 27 to improve some other global performances, such as the consensus of slip ratios. Since K1 is obtained by solving the Riccati equation of a local system, the computation cost is independent of the number of IWM actuators. See Appendix C for the details of the hierarchical LQR algorithm. Notice 11: Similar to Assumption 5, Assumption 10 is just to limit our discussion to homogenous multi-agent systems. Motivated by a recent work of Nguyen and Hara [33], we will develop the H-LQR version for heterogeneous multi-agent systems with physical interaction in our future study.
Augmented model Figure 13. Structure of the hierarchical decentralized LQR using DM-3.

Evaluation of the Proposed Design Models by Carsim/Matlab Co-Simulator
Using the Carsim/Matlab co-simulator introduced in Section 3, this section is to demonstrate the effectiveness of the proposed design models.

Evaluation of Passivity Based Anti-Slip Control Using DM-1
To evaluate the proposed anti-slip control strategy, a test is conducted as follows. As shown in Figure 14a, the EV accelerates from the initial position on the high friction surface (μ = 0.8). It enters the low friction surface (μ = 0.2) at about 5.7 s. The total driving command Tcmd is expressed in Figure  14b. The distribution vector k is updated to minimize the total workload of all wheels [55]. Notice 1: Similar to Assumption 1, Assumption 3 is just to limit our discussion to homogenous multi-agent systems. Motivated by a recent work of Nguyen and Hara [33], we will develop the H-LQR version for heterogeneous multi-agent systems with physical interaction in our future study.

Evaluation of the Proposed Design Models by Carsim/Matlab Co-Simulator
Using the Carsim/Matlab co-simulator introduced in Section 3, this section is to demonstrate the effectiveness of the proposed design models.

Evaluation of Passivity Based Anti-Slip Control Using DM-1
To evaluate the proposed anti-slip control strategy, a test is conducted as follows. As shown in Figure 14a, the EV accelerates from the initial position on the high friction surface (µ = 0.8). It enters the low friction surface (µ = 0.2) at about 5.7 s. The total driving command T cmd is expressed in Figure 14b. The distribution vector k is updated to minimize the total workload of all wheels [55].

Evaluation of GFV Based Stability Analysis Using DM-2
To demonstrate the GFV based stability analysis using DM-2, we conducted wheel velocity control using the system in Figure 9. The EV operates on the high friction surface with μ = 0.85, as Case 1 (Without anti-slip control: Figure 15): The motor torque of each wheel is only distributed from the total driving command, i.e., T i = k i T cmd where k i is the distribution ratio of the ith wheel. Consequently, the slip occurs. Especially, the rear-wheel velocities increase considerably as the vehicle enters the low friction surface.   Case 2 (With passivity-based anti-slip control: Figure 16): The anti-slip control law is implemented based on passivity and the DM-1. By trial-and-error, we select the control gains K a = 100 and K ω = 0.0001. Notice that a small positive value of K ω is enough to assure the strict passivity of the control system. We should not select a big value of K ω which certainly degrades the acceleration performance of the vehicle. In contrast to Case 1, the motor torques are generated properly on different road surfaces. Thus, the wheel velocities are always closed to the vehicle velocity.

Evaluation of GFV Based Stability Analysis Using DM-2
To demonstrate the GFV based stability analysis using DM-2, we conducted wheel velocity control using the system in Figure 9. The EV operates on the high friction surface with μ = 0.85, as shown in Figure 17a. The total driving command Tcmd is expressed in Figure 17b. Similar to the previous sub-section, the distribution vector k is updated to minimize the total workload of all wheels. The driving force observer is designed with the time constant τf = 0.033 [s]. This paper examines the operating points with vxo = 10 (m/s). Considering the high friction surface, it is reasonable to assume that a small slip ratio of 0.05 is maintained constantly when the vehicle velocity changes between vxo − Δv and vxo + Δv where Δv = 2 (m/s). We select the wheel velocity controller with {Kp =52.8; Ki = 528.0}. It is transparent that this controller stabilizes the local transfer function H(s) about the operating The above test cases clarify the merit of the proposed anti-slip strategy using DM-1. Although the physical system of the EV is quite complex and nonlinear, it can be modelled hierarchically with the aggregation and distribution. Based on this model, we can analyze the passivity of the EV dynamics, and design the anti-slip control law that rigorously guarantees the stability of the overall system. Thanks to passivity theory, we design the system with less mathematical and computational effort.
Using the DM-1, the energy storage function of the EV is shown to be S = 1 Given a desired velocity pattern of v x , the energy consumption of the EV can be reduced if the longitudinal velocities of the wheels {rω i } are closed to v x . Therefore, the proposed anti-slip control is not only for improving the safety, but also extending the use of the energy source for the EVs.

Evaluation of GFV Based Stability Analysis Using DM-2
To demonstrate the GFV based stability analysis using DM-2, we conducted wheel velocity control using the system in Figure 9. The EV operates on the high friction surface with µ = 0.85, as shown in Figure 17a. The total driving command T cmd is expressed in Figure 17b. Similar to the previous sub-section, the distribution vector k is updated to minimize the total workload of all wheels. The driving force observer is designed with the time constant τ f = 0.033 [s]. This paper examines the operating points with v xo = 10 (m/s). Considering the high friction surface, it is reasonable to assume that a small slip ratio of 0.05 is maintained constantly when the vehicle velocity changes between v xo − ∆v and v xo + ∆v where ∆v = 2 (m/s). We select the wheel velocity controller with {K p =52.8; K i = 528.0}. It is transparent that this controller stabilizes the local transfer function H(s) about the operating points. Applying the Stability test 1, the stability regions at the operating points of v xo , v xo − ∆v and v xo + ∆v are plotted as the shaded regions in Figure 18a-c, respectively. It is transparent that the points (−4, 0) and (0,0) are always placed in such shaded regions. This means the wheel velocity control system in Figure 9 is stable at the aforementioned operating points.  Figure 18a-c, respectively. It is transparent that the points (−4, 0) and (0,0) are always placed in such shaded regions. This means the wheel velocity control system in Figure 9 is stable at the aforementioned operating points. The above example of stability test unveils a merit of GFV approach using DM-2. Thanks to the rank-1 physical interaction expressed by DM-2 and the GFV theory, the complexity of the stability test is independent of the number of IWM actuators. We can repeat the above test to assure the stability of the control system at different operating points with the vehicle velocity from zero to vxo. This statement is confirmed again by the simulation results summarized in Figure 19. In this simulation, the vehicle is accelerated until 10 s. From 10 s, the total driving command becomes zero, and the vehicle is gradually decelerated due to friction.
In Appendix D, we present several results of the real-time experiments of the wheel velocity control system in Figure 9.      The above example of stability test unveils a merit of GFV approach using DM-2. Thanks to the rank-1 physical interaction expressed by DM-2 and the GFV theory, the complexity of the stability test is independent of the number of IWM actuators. We can repeat the above test to assure the stability of the control system at different operating points with the vehicle velocity from zero to v xo . This statement is confirmed again by the simulation results summarized in Figure 19. In this simulation, the vehicle is accelerated until 10 s. From 10 s, the total driving command becomes zero, and the vehicle is gradually decelerated due to friction.

Evaluation of LQR-Based Slip-Ratio Control Using DM-3
Using the DM-3, we proposed the LQR algorithm for controlling slip-ratio of IWM-EV [17]. This algorithm is verified in the acceleration mode. We also compared the performance of the LQR algorithm with the others, such as sliding mode algorithm and model-following algorithm in [17]. This comparison presented the merits of the proposed LQR: (i) As a hierarchically decentralized algorithm, its computational cost is considerably reduced in comparison with the conventional Figure 19. Simulation results of wheel velocity control using DM-2. In Appendix D, we present several results of the real-time experiments of the wheel velocity control system in Figure 9.

Evaluation of LQR-Based Slip-Ratio Control Using DM-3
Using the DM-3, we proposed the LQR algorithm for controlling slip-ratio of IWM-EV [17]. This algorithm is verified in the acceleration mode. We also compared the performance of the LQR algorithm with the others, such as sliding mode algorithm and model-following algorithm in [17]. This comparison presented the merits of the proposed LQR: (i) As a hierarchically decentralized algorithm, its computational cost is considerably reduced in comparison with the conventional centralized LQR. (ii) Unlike the completely decentralized control algorithms, the physical interaction between the local wheels are properly treated. This gives the degree of freedom to improve the consensus of local wheels' slip ratios. In addition, this can reduce the vibration of slip ratios introduced by the non-flat surfaces. (iii) The optimal control gains can be updated at every control period by the measurements of vehicle motion (vehicle velocity, wheel velocities, driving forces).
For the sake of paper space, we neglect to repeat the comparison between the control approaches in this study. Instead, we only demonstrate the performance of the proposed algorithm in deceleration mode (Figure 20a). A total driving command is given as in Figure 20b. This command is distributed to the local wheels to minimize the total workload. The vehicle runs on high friction surface (µ = 0.8) from the beginning until 12 s. From 12 s, the vehicle starts to decelerate on the low friction surface (µ = 0.2). The proposed LQR algorithm is applied to maintain a slip ratio of λ* = −0.1 during the deceleration period. To this end, we design the LQR algorithm with

Conclusions
This paper suggests that IWM-EV should be treated as a special type of multi-agent-system with physical interaction between the local agents or the locally controlled wheels. Three ways to model To compensate the influence of load transfer between the front and rear wheels during the deceleration period, we select the positive semidefinite matrix Ψ N as follows with ϕ front = ϕ rear = 1.
Simulation results of two test cases are summarized as follows. Case 1 (without slip-ratio control: Figure 21a): The total driving command T cmd in Figure 18b is directly distributed to each wheel. Consequently, the vehicle suffers severe slip during the deceleration periods. The wheel velocities are reduced too fast in comparison with the vehicle velocity. Especially, the wheels reverse their rotation after 14 s.

Conclusions
This paper suggests that IWM-EV should be treated as a special type of multi-agent-system with physical interaction between the local agents or the locally controlled wheels. Three ways to model the IWM-EV are proposed. Each is shown to be suitable for certain control objectives and certain system design approach. From a multi-agent system point of view, the proposed design models, namely DM-1,2,3, contribute some new understandings on IWM-EV dynamics.
Firstly, the DM-1 points out that the IWM-EV naturally has the hierarchical structure: the vehicle body dynamics in the upper-layer, and the local wheel dynamics in the lower-layer. The two layers are connected via the aggregation and distribution channels. Although DM-1 is nonlinear and complex, it is quite helpful for applying passivity theory to clarify the energy storage function of IWM-EV. Passivity theory, therefore, is a possible candidate for improving the safety and optimizing the use of electricity in IWM-EVs.
Secondly, the DM-2 shows another interesting characteristic of IWM-EV. This is the rank-1 physical interaction. Based on DM-2, IWM-EV motion control system can be modelled by a local transfer function H(s) and the interaction matrix. This model is quite useful for applying generalized frequency variable to analyze the stability of the control system. The complexity of stability analysis is independent of the number of IWM actuators.
Last but not least, the DM-3 is a state-space representation of IWM-EV. It is a suitable choice for applying optimal control strategies, such as LQR control of the slip-ratios and driving forces. In addition to attaining the local control objective at each wheel, DM-3 can be utilized to address the additional global control objectives.
This paper demonstrates several control design examples using the proposed DMs. The effectiveness of the DMs and their design approaches are evaluated using Carsim/Matlab cosimulator. For the sake of simplicity, DM-2 and DM-3 are established under the assumptions that the Case 2 (with LQR-based slip-ratio control: Figure 21b): The motor torque of each wheel is the summary of two signals: the command distributed from T cmd expressed in Figure 18b; and the output of the LQR-based slip-ratio controller. Thanks to the proposed method, the wheel velocities are slightly smaller than the vehicle velocity during the deceleration period. This means the safety and controllability of the vehicle are successfully maintained. In addition, the efficiency of electric energy utilization is also improved by maintaining a small slip ratio.

Conclusions
This paper suggests that IWM-EV should be treated as a special type of multi-agent-system with physical interaction between the local agents or the locally controlled wheels. Three ways to model the IWM-EV are proposed. Each is shown to be suitable for certain control objectives and certain system design approach. From a multi-agent system point of view, the proposed design models, namely DM-1,2,3, contribute some new understandings on IWM-EV dynamics.
Firstly, the DM-1 points out that the IWM-EV naturally has the hierarchical structure: the vehicle body dynamics in the upper-layer, and the local wheel dynamics in the lower-layer. The two layers are connected via the aggregation and distribution channels. Although DM-1 is nonlinear and complex, it is quite helpful for applying passivity theory to clarify the energy storage function of IWM-EV. Passivity theory, therefore, is a possible candidate for improving the safety and optimizing the use of electricity in IWM-EVs.
Secondly, the DM-2 shows another interesting characteristic of IWM-EV. This is the rank-1 physical interaction. Based on DM-2, IWM-EV motion control system can be modelled by a local transfer function H(s) and the interaction matrix. This model is quite useful for applying generalized frequency variable to analyze the stability of the control system. The complexity of stability analysis is independent of the number of IWM actuators.
Last but not least, the DM-3 is a state-space representation of IWM-EV. It is a suitable choice for applying optimal control strategies, such as LQR control of the slip-ratios and driving forces. In addition to attaining the local control objective at each wheel, DM-3 can be utilized to address the additional global control objectives.
This paper demonstrates several control design examples using the proposed DMs. The effectiveness of the DMs and their design approaches are evaluated using Carsim/Matlab co-simulator. For the sake of simplicity, DM-2 and DM-3 are established under the assumptions that the local agents are homogeneous (Assumptions 1 and 3). To overcome this real limitation, in future study we will employ the multi-agent theories to deal with the heterogeneous of the local wheels.
This paper only focuses on the dynamical models for vehicle motion control. Therefore, we do not consider the models of other parts of the EV systems, such as the motor model, the inverter model, and the energy source model. Motivated by several recent works on battery modelling [37][38][39][40], we will also develop the full model of the EVs using Matlab-Carsim co-simulator. Such simulation models will be utilized to integrate motion control, motor control, and optimal energy management for IWM-EVs. Funding: This study is supported in part by the Toyota Technological Institute (Aichi, Japan).

Acknowledgments:
We would like to thank Shinji Hara (Tokyo Institute of Technology) who gave kind advice in the theoretical background of this paper. We also thank Hori-Fujimoto Lab (The University of Tokyo) for the experimental electric vehicle presented in Appendix A.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1 summarizes some main parameters of the simulation model which is based on the electric vehicle prototype shown in Figure A1.  local agents are homogeneous (Assumptions 5 and 10). To overcome this real limitation, in future study we will employ the multi-agent theories to deal with the heterogeneous of the local wheels. This paper only focuses on the dynamical models for vehicle motion control. Therefore, we do not consider the models of other parts of the EV systems, such as the motor model, the inverter model, and the energy source model. Motivated by several recent works on battery modelling [37][38][39][40], we will also develop the full model of the EVs using Matlab-Carsim co-simulator. Such simulation models will be utilized to integrate motion control, motor control, and optimal energy management for IWM-EVs.

Acknowledgments:
We would like to thank Shinji Hara (Tokyo Institute of Technology) who gave kind advice in the theoretical background of this paper. We also thank Hori-Fujimoto Lab (The University of Tokyo) for the experimental electric vehicle presented in Appendix D.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1 summarizes some main parameters of the simulation model which is based on the electric vehicle prototype shown in Figure A1.

Appendix B
The GFV φ + + + + = + + + To get an optimal feedback gain of the form (44), the weighting matrices are selected as where the matrices (Q 1 , R 1 , H 1 ) are set for the lower layer, (Q g1 , R g1 ) and (Q g2 , R g2 ) are set for the upper layer. The control gains are obtained by the three-step procedure as follows Step 1 (Lower-layer): Select the weighting matrices (Q 1 , R 1 , S 1 ) such that Q 1 and S 1 are positive semidefinite, and R 1 is strictly positive definite. Then, solve the differential Riccati equation with the boundary condition P 1 t f = P 1 f = H 1 .
Step 2 (Upper-layer): Select the positive semidefinite matrix Ψ N and the positive definite matrices R g1 and R g2 to obtain the weighting matrix R as in (45). Then, the two matrices Q g1 and Q g2 are set as following such that Q g1 is positive semidefinite: Step 3 (Feedback gain): Finally, the feedback gains are calculated as The above procedure is a kind of partial inverse LQR method that provides the hierarchical control structure in Figure 13, since the choices of Q g1 and Q g2 are not completely free but restricted as in (45). The optimality of the feedback gains is shown by Proposition 1 in [17].

Appendix D
In this paper, we present the general design models for the electric vehicle driven by N IWM actuators. All the simulations are performed by using a 4-wheel-model. We actually implemented several motion control approaches to a real electric vehicle, such as lateral motion control in [3], and longitudinal motion control in [17]. It is the electric vehicle COMS shown in Figure A2. It has four wheels but only two rear wheels are driven by IWMs of IPMSM type with the maximum power of 2 kW. The motor drives and inverters are provided by Myway. The vehicle is equipped with RTK-GPS receiver produced by Hemisphere, the gyroscopes, and the accelerometers. Through the fusion of GPS receiver, gyroscopes, and accelerometers, the vehicle velocity can be obtained. The encoders are installed at the wheels to measure the rotational speed. The heart of the system is a RT-Linux computer which processes the control algorithm and stores the experimental data.
In Figure A3, we demonstrate an experimental result of the wheel velocity control system in Figure 9. This was an autonomous driving test in which the vehicle is required to accelerate to the velocity of 5 m/s and then, maintaining this velocity constantly. To this end, the driving command is given by the upper-layer motion control layer. To design the PI controller for controlling the wheel velocity, we apply the Stability condition 8 in case N = 2. Thanks to this condition, we can stabilize the overall system and smoothly attain both the global motion and local motions of the vehicle.