A Steering-Following Dynamic Model with Driver’s NMS Characteristic for Human-Vehicle Shared Control

Abstract

For investigating driver characteristic as well as control authority allocation during the process of human–vehicle shared control (HVSC) for an autonomous vehicle (AV), a HVSC dynamic mode with a driver’s neuromuscular (NMS) state parameters was proposed in this paper. It takes into account the driver’s NMS characteristics such as stretch reflection and reflex stiffness. By designing a model predictive control (MPC) controller, the vehicle’s state feedback and driver’s state are incorporated to construct the HVSC dynamic model. For the validation of the model, a field experiment was conducted. The vehicle state signals are collected by V-BOX, and the driver’s state signals are obtained with the electromyography instrument. Subsequently, the hierarchical least square (HLS) parameter identification algorithm was implemented to identify the parameters of the model based on the experimental results. Moreover, the Unscented Kalman Filter (UKF) was utilized to estimate the important NMS parameters which cannot be measured directly. The experimental results showed that the model we proposed has excellent accuracy in characterizing the vehicle’s dynamic state and estimating the driver’s NMS parameter. This paper will serve as a theoretical basis for the new control strategy allocation between human and vehicle for L3 class AVs.

1. Introduction

Significant progress has been conducted in the field of perception, decision making, path planning and control authority for autonomous vehicles (AVs) in the past few decades. However, it is still far away from widespread market penetration because of its unmatured technology. As a promising technique aimed at improving the safety of AV, the concept of “human–vehicle shared control” (HVSC) has achieved a great deal of research effort in recent years [1]. For HVSC, the automation and driver take control authority of the vehicle sequentially after each other, or co-pilot the vehicle simultaneously [2].

The dynamic model is of great significance for the understanding of the intrinsic and developing control authority of HVSC. In this way, many existing scholarly works have demonstrated the contributions associated with the dynamic model in the HVSC domain. For example, the biomechanical properties of the driver for a vehicle steering task is considered a novel technique of estimating the admittance and total response [3]. By investigating the influence of haptic aids on the pilot’s NMS response, an online estimator of the time-varying NMS dynamic method based on Recursive Least Squares is proposed [4,5]. On the basis of the electric power-assisted steering system, Marouf A et al. [6] proposed a time-varying method based on a sliding mode observer and the exponentially weighted recursive least squares (EWRLS) algorithm for the estimation of the viscoelastic properties in the driver’s arm. The characterization of driver NMS dynamics, implemented to analyze the influence of active and passive steering tasks and hand positions, is also highly sought after [7]. Subsequently, in order to further reflect the relationship between driver and vehicle, the novel Queuing Network is integrated into the driver’s NMS dynamic [8]. Similarly, by reducing the internal cost function of the driver, a human–vehicle steering control strategy with the driver’s NMS dynamic is proposed [9]. In this way, an intra-driver variability human-like steering model inspired the research boom [10]. Soon after, Yanbo Z et al. [11] proposed a shared steering controller by imitating the driver’s characteristics of the NMS system. Furthermore, other physiological characteristics including the driver’s vision and blood volume pulse, which are regarded as vital indexes, are utilized to demonstrate the mechanism of driving behavior and construct the human–vehicle shared control model [12,13].

When the driver executes the operation of following the leading car, the driver’s leg model plays a vital role in HVSC. In regard to this field, the pedal drive unit is applied to the human–machine interface analysis experiment by investigating the physiological characteristics of the driver’s leg [14]. It is noteworthy that the mechanism of the individual muscle-tendon biomechanics metabolism also is parsed, and then a Hill-based human leg model with a muscle and joint is proposed [15]. In further research, a put-forward two-dimensional leg model with a knee, along with ankle passive and quasi-passive elements is optimized [16]. By investigating the influence of muscle force during flexion–extension movement, a Hill-based electromyography (EMG)-driven model of the elbow is proposed [17].

In light of our literature review, there are only a few studies involved in the following aspects. The current driver dynamic models with NMS characteristic are only represented by a single steering arm. The leg NMS dynamic of the driver was not taken into account thoroughly. Particularly in the condition of a roundabout or overtaking, the driver executes the sequential steering–following operation. Therefore, the relative influence between arm and leg NMS characteristic should not be neglected for practical application of the dynamic model. For L3 autonomous vehicle (AV), the driver is required to take over the control authority of the vehicle when the automation system fails to deal with the emergency situation. Most of time, the driver will be disengaged from the human–vehicle–road closed loop. Therefore, frequently, discrete interactions between the driver and vehicle will cause the degradation of the human’s cognitive ability and driving skills. Consequently, it will affect the driver’s NMS dynamic properties, such as viscoelastic resistance coefficient and movement feedback coefficient. Moreover, the driver’s NMS model was not analyzed by employing the EMG signals of the driver for validation.

In this paper, a HVSC dynamic model incorporating the arm and leg NMS characteristics of the driver is proposed. The EMG signals are collected and then muscle force is decomposed based on the Hill muscle contraction theory. The hierarchical least square method is applied to identify the off-line parameters of the dynamic model. Thereafter, the Unscented Kalman Filter (UKF) observer is designed to estimate the NMS characteristic parameters. The structure of the article is organized as follows: Section 2 illustrates the schematic of HVSC model. Section 3 describes the scenario and the procedure of the field experiment. The offline identification of the model’s parameter and online estimation by means of UKF observer are demonstrated in Section 4. The proposed dynamic model is validated using experimental results in Section 5. A brief summary and discussion are present in the last section.

2. Materials and Methods

By designing a model predictive control (MPC) controller with vehicle’s state feedback, a HVSC dynamic model is constructed in this paper. Specifically, the function of the MPC controller is to keep tracking the trajectory to generate the anticipated steering wheel angle ${\widehat{\delta }}_{sw}$ and gas pedal position ${\widehat{\theta }}_{ped}$. After receiving the anticipated angle ${\widehat{\delta }}_{sw}$, the driver’s NMS model produces the actual steering angle ${\delta }_{sw}$ and gas pedal position ${\theta }_{ped(f)}$ signal feedback to controller backward. The NMS dynamic model of the driver is composed of 11 components. The schematic of the HVSC dynamic model is illustrated as Figure 1.

2.1. Subsection

The dynamic response of the driver’s arms applied at the steering wheel can be represented by a parallel spring and damper, as shown in

Table 1. Parameters of driver neuromuscular (NMS) model.

Symbol	Description	Value
Jst	Steering wheel inertia	0.1395 (kg m2)
Bst	Steering wheel damper	1.56 (Nm s/rad)
Kst	Steering wheel spring	2.29 (Nm/rad)
Ka	Active stiffness	−1920 (Nm/rad)
Jdr	Arm inertia	0.172 (kg m2)
Bdr	Arm damping	1.032 (Nm s/rad)
Kdr	Arm spring	60.707 (Nm/rad)
Kr	Arm reflection stiffness	3.423 (Nm s/rad)
Br	Arm reflection damping	1.69 (Nm/rad)
τdel	Inherent time delay	0.04 (s)
Kpos	Stretch feedback gain	0.52 (Nm/rad)
Ktend	Tendon stiffness	2799 (Nm/rad)
f0	Eigen-frequency	1.1 (Hz)
Iseg	Motion of the inertia	105.5 (g/m2)
Kf	GTO feedback gain	1.18 (Nm/rad)
Kvel	Stretch velocity feedback gain	40.4 (Nm/rad)
Kint	Co-contraction stiffness	334 (Nm/rad)
Bint	Co-contraction damping	19.4 (Nm/rad)
Kcon	Contact elasticity	1033 (Nm/rad)
Bcon	Contact viscosity	11.2 (Nm /rad)

. The passive damping and stiffness of the arms represent the properties of the muscles, skin and tendons in relaxation, which can be simplified as arm inertia $J_{dr}$, arm damping $B_{dr}$ and arm stiffness $K_{dr}$ systems [18]. The steering system, including the steering wheel, column, rack and pinion, uprights and wheels and tires can also be represented as simple inertia $J_{st}$, damper $B_{st}$ and spring $K_{st}$ systems. In addition, the torque feedback term $K_a$ is inferred to the torque $M_t$ obtained from the slip angle at the front axle [19].

The steering mechanism dynamics and steering arm dynamics are coupled at the point where the driver grips the steering wheel. The equation of motion and corresponding transfer function is obtained as follows:

$$(J_{dr}+J_{st}){\ddot{\delta }}_{\mathrm{sw}}+(B_{dr}+B_{st}){\dot{\delta }}_{\mathrm{sw}}+(K_{dr}+K_{st}){\delta }_{\mathrm{sw}}=T_m-\frac{M_T}{n_{rsw}}$$

(1)

$$G_s(s)=1/[(J_{dr}+J_{st})s^2+(B_{dr}+B_{st})s+(K_{dr}+K_{st})]$$

(2)

where $T_m$ represents the total muscle torque applied at the steering wheel; $n_{rsw}$ and ${\delta }_{\mathrm{sw}}$ indicate the steering gear ratio and wheel steering angle respectively.

The spindle feedback muscle information is transferred to the alpha-motor neurons to maintain muscle length. The reflex control of muscle length can cause a change in the intrinsic passive stiffness and damping of the muscle. The value of the reflex gain varies depending on muscle activation [20]. The reflex control dynamic can be expressed as Equation (3):

$$H_r(s)=[{\omega }_c(s{B}_r+K_r)e^{-s{\tau }_{del}}]/(s+{\omega }_c)$$

(3)

where ${\omega }_c$ refers to the cutoff frequency; $B_r$ indicates the reflection damping.

Noteworthily, the inherent time delay can be approximated as a first-order transfer function:

$$e^{-s{\tau }_{del}}\approx 1/(1+{\tau }_{del}s)$$

(4)

The steering torque command is generated by a reference model, which is, essentially, a learned inverse model of the arms, steering and vehicle. The following reference model dynamic is extremely appropriate for real driver behavior [20].

$$\widehat{G}(s)=\frac{0.067}{0.028s^2+0.083s+1}$$

(5)

2.2. Leg NMS Dynamic Model

The admittance $H_{adm}$ can be parameterized by a linear NMS model which represents the dynamics of the ankle–foot interaction with the gas pedal in Figure 1. The inputs of the NMS model are the disturbing torque $T_{dis}$ and the anticipated gas pedal position ${\widehat{\theta }}_{ped}$. The measurable model outputs are $T_c$ and ${\theta }_c$ [21]. The model structure is an extension of NMS models proposed in previous research [22].

The inertia component can be described by the endpoint inertia of the limb $I_{seg}$ in the form of a transfer function as follows:

$$H_{seg}(s)=\frac{1}{I_{seg}{s}^2}$$

(6)

The muscle contact dynamic is described as:

$$H_c(s)=k_{con}+b_{con}s$$

(7)

The intrinsic component is defined by the activated muscle stiffness $k_{\mathrm{int}}$ and damping $b_{\mathrm{int}}$ in the form of transfer function:

$$H_{\mathrm{int}}(s)=k_{\mathrm{int}}+b_{\mathrm{int}}s$$

(8)

The reflexive component can be deduced as follows:

$$H_{act}(s)=\frac{1}{(\frac{1}{{\omega }_0^2})s^2+\frac{2\beta }{{\omega }_0}s+1}$$

(9)

The muscle golgi tendon organs (GTO) dynamic can be demonstrated as:

$$H_{got}(s)=k_f{e}^{-{\tau }_{del}s}$$

(10)

The muscle spindle dynamic is expressed as:

$$H_{ms}(s)=(K_{pos}+K_{vel}s)e^{-s·{\tau }_{del}}$$

(11)

where the $K_{pos}$ and $K_{vel}$ represent the gains of the monosynaptic stretch and stretch velocity feedback respectively.

As a serial elastic element to the muscles, tendons dynamic can be expressed as:

$$H_{tend}(s)=\frac{1}{-k_{tend}}$$

(12)

where $k_{tend}$ represents the stiffness of the tendon.

Noteworthily, the motion of the inertia $I_{seg}$ is the result of the sum of $T_c(s)$ and $T_{mus}(s)$. The pedal–torque can be combined with the muscles, partially. Then, the motion of the limb is deduced as:

$${\theta }_{\lim b}(s)=H_{seg}(s)[T_c(s)+T_{mus}(s)]$$

(13)

Pedal rotations do not only cause joint rotations but also small displacements of the skin or soft tissue. This effect is described by contact dynamics as follows:

$$T_c(s)=H_c(s)[{\theta }_c(s)-{\theta }_{\lim b}(s)]$$

(14)

2.3. Vehicle Dynamic Model

During the process of HVSC, the lateral and yaw motion of the vehicle are critical issues. A simplified linear two-degree-of-freedom (2DOF) lateral dynamic model is widely utilized to investigate the characteristic of human–automation interaction [23]. The geometric relationship and parameters of the vehicle dynamic model are listed in Figure 2 and

Table 2. Parameters of the vehicle dynamic model.

Symbol	Description	Value
Iz	Yaw inertia	4192 (kg m2)
vx	Longitudinal speed	2.29 (m/s)
k1	Cornering stiffness of the front tires	93360 (N/rad)
k2	Cornering stiffness of the rear tires	57340 (N/rad)
a	Front axle length	1.18 (m)
b	Rear axle length	1.58 (m)
m	Body mass	1498 (kg)
xl	Look-ahead distance	12 (m)
nrsw	Steering gear ratio	18

:

The equations of the lateral and yaw moment are denoted as:

$$\{\begin{array}{l}{\dot{v}}_y=-\frac{(k_1+k_2)}{m{v}_x}{v}_y-(v_x-\frac{(b{c}_r-a{c}_f)}{m{v}_x}){\omega }_r+\frac{k_1}{m}{\delta }_{sw}\hfill \\ {\dot{\omega }}_r=\frac{(b{k}_2-a{k}_1)}{I_z{v}_x}{v}_y-\frac{(b^2{c}_r+a^2{c}_f)}{I_{\mathrm{z}}{v}_x}{\omega }_r+\frac{a{k}_1}{I_z}{\delta }_{sw}\hfill \end{array}$$

(15)

where v_y represents lateral speed; ω_r infers to yaw rate. The lateral offset e_y is described as:

$${\dot{e}}_y=v_y+x_l{\omega }_r+v_x\psi$$

(16)

where $\psi$ represents the yaw angle.

2.4. MPC Controller

In order to incorporate the vehicle dynamics with road trajectory, a great deal of investigation of trajectory tracking was conducted. Inspired by appropriate literature, the MPC controller is implemented in this section for strong robustness and stability in trajectory tracking [24]. The main steps of the MPC controller are given in Algorithm 1:

Algorithm 1. MPC controller.

Input_U: {Vehicle state variable $V_{cur}$; Road trajectory} Step1: Construct vehicle predictive model: $\xi (k|t)=\left[\begin{array}{c}V(k|t)\\ u(k-1|t)\end{array}\right]$ $\{\begin{array}{l}\xi (k+1|t)={\mathrm{A}}_{k,t}\xi (k|t)+{\mathrm{B}}_{k,t}\Delta U(k|t)\hfill \\ V(k|t)=C_{k,t}\xi (k|t)\hfill \end{array}$ $y(t)={\mathsf{\Psi }}_t\xi (t|t)+{\mathsf{\Theta }}_t\Delta U(t)$ Step2: Calculate cost function: $\begin{array}{cc}\hfill J(t)& ={\displaystyle \sum _{i=1}^{N_p}{\Vert V(t+i|t)-V_{ref}(t+i|t)\Vert }_Q^2}\hfill \\ & +{\displaystyle \sum _{i=1}^{N_{c-1}}{\Vert △u(t+i|t)\Vert }_R^2+\rho {\epsilon }^2}\hfill \end{array}$ Step3: Set constraints: $\Delta {\mathrm{U}}_{\min }\le \Delta {\mathrm{U}}_t\le \Delta {\mathrm{U}}_{\max }$ ${\mathrm{U}}_{\min }\le \mathrm{A}\Delta {\mathrm{U}}_t+{\mathrm{U}}_t\le {\mathrm{U}}_{\max }$ $y_{hc,\min }\le y_{hc}\le y_{hc,\max }$ $y_{sc,\min }-\epsilon \le y_{sc}\le y_{sc,\max }+\epsilon$ $\epsilon >0$ Output_y: {Anticipated steering wheel angle ${\widehat{\delta }}_{sw}$; Anticipated pedal position ${\widehat{\theta }}_{ped}$}

2.5. Model Integration

Considering that the vehicle dynamic model is constructed in the form of state-space equation, the driver’s NMS dynamic expressed as transfer functions should be converted to the state equation. From the aspect of theoretic analysis, the n-order transfer function can be decomposed into the product of several one-order transfer functions. Then the system state variable diagram is formed in series by one-order transfer functions [25]. In this way, the HVSC state equation form the block diagram can be constructed as shown in Figure 3.

A total of 17 state variables are defined in Figure 3.The physical sense of ${\delta }_{\mathrm{sw}}$, ${\dot{\delta }}_{sw}$, $T_{aca}$, $T_{\xi c}$, $T_{rc}$, $T_{cca}$, ${\theta }_{\lim b}$, ${{\theta }^{\prime }}_{\lim b}$, $T_{ccl}$, $T_{\mathrm{int}}$, $T_{acl}$, ${T^{\prime }}_{acl}$, $T_{got}$ and $T_{sp}$ represents the steering wheel angle, steering wheel angular rate, arm active contraction torque, reference contraction torque, contact contraction torque, reflex contraction torque, leg motion, leg motion rate, leg contact contraction torque, leg intrinsic contraction torque, leg active contraction torque, leg active reference torque, leg GTO feedback torque and leg spindle feedback torque respectively. Noteworthily, the physical value of some state variables cannot be obtained directly from the sensor. Therefore, a state observer based on UKF is proposed to estimate these variables. The relationship between state variables and measurable variables was elaborated on in the following content [26].

First, the muscle contraction force $F_i$ is deduced based on Hill muscle theory as follows:

$$\begin{array}{cc}\hfill F_i& =F_{\mathrm{ce}}+F_{\mathrm{pe}}=F_{\max }[a(i)f_a(\epsilon )+f_p(\epsilon )]\hfill \\ & =F_{\max }[a(i)e^{-{[2.7272778\times \ln (\epsilon +1)]}^2}+0.03(e^{7\epsilon }-1)]\hfill \end{array}$$

(17)

where F_ce represents muscle contraction force; F_pe infers to the elastic force; $F_{\max }$ indicates the muscle maximum contraction force; $a(i)$ is the muscle activation; $f_a(\epsilon )$ and $f_p(\epsilon )$ infer to the active and passive contraction of the muscle; $\epsilon$ represents deformation degree of the muscle during contraction.

Afterward, the weighted analysis method is used to calculate the arm active contraction torque $T_{aca}$ and leg active contraction torque $T_{acl}$ applied to the steering wheel and pedal respectively as follows:

$$\{\begin{array}{l}{T}_{aca}=(F_1{n}_1+F_2{n}_2+\cdots +F_i{n}_i)/R_a;\stackrel{}{}i=6\\ T_{acl}=(F_1{n}_1+F_2{n}_2+\cdots +F_i{n}_i)/R_l;\stackrel{}{}i=4\end{array}$$

(18)

where $n_i$ infers to the weight coefficient of muscle force; R_a and R_l represent the radius of the steering wheel and pedal rotation.

For the arm dynamic model, the arm reflex contraction torque $T_{rc}$ can be deduced based on $T_{aca}$ and assist torque $T_{as}$ applied on the steering wheel, which can be express as Equation (19):

$$T_{rc}=T_{aca}-T_{as}$$

(19)

The arm reference contraction torque $T_{\xi c}$ and arm contact contraction torque $T_{cca}$ can be described as follows:

$$T_{\xi c}+T_{cca}+K_a({\theta }_{sw}-{\widehat{\theta }}_{sw})=T_{dr}$$

(20)

$$T_{cca}=T_{sw}-T_{dis}$$

(21)

where $K_a$ represents the active stiffness; $T_{dr}$ infers to the steering wheel torque; $T_{sw}$ infers to the torque applied on the steering wheel; $T_{dr}$ indicates the steering torque of vehicle; $T_{dis}$ is the torque perturbation from the sensor.

As shown in Figure 3, the feedback torque $T_{feed}$ can be calculated as steering angle feedback and acquired steering angle, denoted as Equation (22)

$$T_{feed}={\delta }_{\mathrm{sw}}{K}_e+\stackrel{·}{{\delta }_{\mathrm{sw}}}{B}_e+G_a/n_{rsw}$$

(22)

where $T_{feed}$ is the steering wheel torque feedback; $K_e$ is the gain coefficient of steering angle; $B_e$ represents the gain coefficient of steering speed; $G_e$ represents the kingpin caster angle gain coefficient.

For the leg dynamic model, the leg contact contraction torque $T_{ccl}$ can be expressed by the pedal torque $T_{ped}$:

$$T_{ccl}=T_{ped}-T_{dis}$$

(23)

The leg intrinsic contraction torque $T_{\mathrm{int}}$ reflects the internal characteristics of leg muscle, which can be represented according to the relationship in Figure 3.

$$T_{\mathrm{int}}=T_{acl}+K_{tend}{\theta }_{tend}$$

(24)

The GTO feedback torque $T_{got}$ the spindle feedback torque and $T_{sp}$ can be deduced as follows:

$${\theta }_{\lim b}={\widehat{\theta }}_{ped}-T_{got}\tau -T_{sp}\tau$$

(25)

where ${\widehat{\theta }}_{ped}$ represent the anticipate gas pedal position.

As described above, the driver’s NMS model expressed as a stated equation can be combined with the vehicle dynamic model expressed in Equation (15). In this way, the final HVSC dynamic model can be demonstrated as an augmented matrix.

$$\{\begin{array}{l}\stackrel{·}{x(t)}=\left[\begin{array}{cc}{A}_1& 0\\ 0& A_2\end{array}\right]x(t)+\left[\begin{array}{c}{B}_1\\ B_2\end{array}\right]u(t)\hfill \\ y(t)=\left[\begin{array}{cc}{C}_1& 0\\ 0& C_2\end{array}\right]x(t)\hfill \end{array}$$

(26)

$$A_1=\left[\begin{array}{cccccccccc}0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ a_{21}& a_{22}& a_{23}& a_{24}& 0& 0& 0& 0& 0& a_{110}\\ a_{31}& 0& a_{33}& 0& 0& a_{36}& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& a_{54}& a_{55}& 0& 0& 0& 0& 0\\ a_{61}& 0& 0& 0& 0& a_{66}& 0& 0& 0& 0\\ a_{71}& 0& 0& 0& 0& 0& a_{77}& a_{78}& 0& 0\\ a_{81}& 0& 0& 0& 0& 0& a_{87}& a_{88}& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& a_{98}& 0& a_{910}\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\end{array}\right]$$

(27)

$$A_2=\left[\begin{array}{cccccccc}0& 1& 0& 0& 0& 0& 0& 0\\ a_{1211}& 0& a_{1213}& a_{1214}& 0& 0& 0& 0\\ a_{1311}& 0& a_{1313}& 0& 0& 0& 0& 0\\ a_{1411}& 0& 0& a_{1414}& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& a_{1614}& a_{1615}& a_{1616}& a_{1617}& a_{1618}\\ a_{1711}& 0& 0& a_{1714}& 0& 0& a_{1717}& 0\\ 0& 0& 0& a_{1814}& 0& 0& 0& a_{1818}\end{array}\right]$$

(28)

$$B_1={\left[\begin{array}{cccccccccc}0& n_2& n_3& 0& n_5& n_6& 0& 0& 0& 0\end{array}\right]}^T$$

(29)

$$B_2={\left[\begin{array}{cccccccc}0& 0& 0& 0& 0& n_{16}& 0& 0\end{array}\right]}^T$$

(30)

$$C_1=\left[\begin{array}{c}\begin{array}{cccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ K_a& 0& 1& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\end{array}\right]$$

(31)

$$C_2=\left[\begin{array}{c}\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& m_{43}& m_{44}& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& -\tau & 0\\ 1& 0& 0& 0& 0& 0& 0& -\tau \end{array}\end{array}\right]$$

(32)

$$D={\left[\begin{array}{cccccccc}0& 0& 0& \cdots \cdots & 0& 0& 0& 0\end{array}\right]}_{1*16}^T$$

(33)

where $a_{{}_{21}}=-(\frac{K_{dr}-K_{st}+K_a-k_{fd}{K}_{\mathrm{e}}+1}{J_{dr}+J_{st}})$; $a_{22}=-(\frac{B_{dr}+B_{st}-B_e{k}_{fd}}{J_{dr}+J_{st}})$; $a_{{}_{23}}=\frac{1}{J_{dr}+J_{st}}$; $a_{{}_{24}}=\frac{1}{J_{{}_{dr}}+J_{{}_{st}}}$; $a_{210}=-\frac{G{}_{a}k{}_{fd}}{(J_{dr}+J_{st})n_{rsw}}$; $a_{31}=-\frac{W_c{B}_r}{I_r}$; $a_{33}=1/I_r$; $a_{36}=1/I_r$; $a_{54}=-1/J_k$; $a_{55}=-B_k/J_k$; $a_{61}=-(W_c{K}_r-W_c{}^2{B}_r)$; $a_{66}=W_{\mathrm{c}}$; $n_2=\frac{W_c{B}_r+K_a}{J_{dr}+J_{st}}$; $n_3=\frac{W_c{B}_r}{I_r}$; $n_5=M_{\mathrm{k}}$; $n_6=W_c{K}_r-W_c{}^2{B}_r$; $a_{71}=\frac{k_1}{m}$; $a_{77}=-\frac{(k_1+k_2)}{m{v}_x}$; $a_{78}=v_x-\frac{(b{c}_r-a{c}_f)}{m{v}_x}$; $a_{81}=\frac{a{K}_1}{I_z}$; $a_{87}=\frac{(b{k}_2-a{k}_1)}{I_z{v}_x}$; $a_{88}=-\frac{(b^2{c}_r+a^2{c}_f)}{I_{\mathrm{z}}{v}_x}$; $a_{98}=x_l$; $a_{910}=v_x$; $a_{1211}=\frac{I_{ped}-K_{tend}}{I_{seg}}$; $a_{1213}=\frac{I_{ped}}{I_{seg}{b}_{con}}$; $a_{1214}=\frac{K_{tend}}{I_{seg}{b}_{\mathrm{int}}}$; $a_{1311}=I_{ped}$; $a_{1313}=\frac{I_{ped}-K_{con}}{b_{con}}$; $a_{1411}=K_{tend}$; $a_{1414}=-\frac{K_{\mathrm{int}}+K_{tend}}{b_{\mathrm{int}}}$; $a_{1614}=-\frac{a_1{\omega }_o^2}{b_{\mathrm{int}}}$; $a_{1615}={\omega }_o^2$; $a_{1616}=-2\beta {\omega }_o$; $a_{1617}=-{\omega }_o^2\tau$; $a_{1618}=-{\omega }_o^2\tau$; $a_{1711}=-K_f{K}_{tend}$; $a_{1714}=\frac{K_f{K}_{tend}}{b_{\mathrm{int}}}$; $a_{1717}=-\tau$; $a_{1814}=\frac{K_{pos}-K_{vel}}{\tau b_{\mathrm{int}}}$; $a_{1818}=-\frac{1}{\tau }$; $m_{43}=-K_{tend}$; $m_{44}=-K_{tend}$. $x^T={[\begin{array}{cc}\begin{array}{lllllllll}{\delta }_{sw}\hfill & {\dot{\delta }}_{sw}\hfill & T_{aca}\hfill & T_{\xi c}\hfill & T_{rc}\hfill & T_{cca}\hfill & {\theta }_{\lim b}\hfill & {{\theta }^{\prime }}_{\lim b}\hfill & T_{ccl}\hfill \end{array}& \begin{array}{lllllllll}{T}_{\mathrm{int}}\hfill & T_{acl}\hfill & {T^{\prime }}_{acl}\hfill & T_{got}\hfill & T_{sp}\hfill & V_y\hfill & {\omega }_{\gamma }\hfill & e_y\hfill & \psi \hfill \end{array}\end{array}]}_{}^T$; $y^T=[\begin{array}{cc}\begin{array}{llllllllll}{\delta }_{sw}\hfill & {\dot{\delta }}_{sw}\hfill & T_{aca}\hfill & T_{dr}\hfill & 0\hfill & T_{sw}\hfill & V_y\hfill & {\omega }_{\gamma }\hfill & e_y\hfill & \phi \hfill \end{array}& \begin{array}{llllllll}{\theta }_{\lim b}\hfill & {\theta }_{\lim b}^{\prime }\hfill & T_{ped}\hfill & {\theta }_{tend}\hfill & T_{acl}\hfill & 0\hfill & 0\hfill & {\theta }_{ped}\hfill \end{array}\end{array}]$; $\mathrm{U}=[\begin{array}{cc}{\delta }_{sw}& {\theta }_{ref}\end{array}]$.

3. Model Parameter Identification and State Observer

Substantially, the HVSC dynamic model established in Section 2 is a multi-input multi-output (MIMO) system with continuous time-varying parameters. The research on the identification algorithm for the linear time-varying system was addressed in numerous works [27,28]. The hierarchical least square (HLS) parameter identification algorithm, which has stronger convergence than recursive least squares algorithm, is utilized in this paper. The principal process of the identification can be summarized as follows.

3.1. Predicting Process

(a)

Initializing the dynamic model: the original HVSC dynamic model demonstrated as continuous state variables should be discretized:

$$\{\begin{array}{l}x(t+1)=Ax(t)+Bu(t)\hfill \\ y(t)=Cx(t)+Du(t)+v(t)\hfill \end{array}$$

(34)

where $v(t)$ is white noise vector with zero mean value; A and B indicate the parameter vector identified in Equation (26).

(b)

Converting the transfer function with the Laplace transform: The relation between input and output is expressed as follows:

$$y(s)=\frac{\lambda (s)}{\eta (z)}u(t)+v(t)=s(t)+v(t)$$

(35)

where $\lambda (z)$ and $\eta (z)$ indicate the molecular and denominator polynomial respectively:

(c)

Recombining the equation with hierarchical iteration:

$$\{\begin{array}{l}\varphi (t)={[-u(t-1),-u(t-2),\cdots ,-u(t-10)]}^T\\ \epsilon (t)=[-s(t-1),-s(t-2),\cdots ,-s(t-10)]\end{array}$$

(36)

Using these definitions, Equation (28) can be deduced as:

$$y(t)=\epsilon (t)\lambda +{\eta }^T\varphi (t)+v(t)$$

(37)

3.2. Identifying Process

(a)

Calculating quadratic cost functions:

$$J_i(\lambda )=\underset{j=1}{\overset{10}{{\displaystyle \mathsf{\Sigma }}}}\parallel y(j)-\epsilon (j)\lambda -{\eta }^T\varphi (t){\parallel }^2$$

(38)

(b)

Identifying the parameter vector $\lambda$ and parameter matrix $\eta$:

$$\{\begin{array}{l}\widehat{\lambda }(t)=\widehat{\lambda }(t-1)+L_1(t)[y(t)-{\widehat{\eta }}^{\mathrm{T}}(t-1)\varphi (t)-\epsilon (t)\widehat{\lambda }(t-1)]\hfill \\ \widehat{\eta }(t)=\widehat{\eta }(t-1)+L_2(t){[y(t)-{\widehat{\eta }}^T(t-1)\epsilon (t)-{\epsilon }^T(t)\widehat{\eta }(t-1)]}^T\hfill \end{array}$$

(39)

$$L_i(t)=P_i(t){\epsilon }^T(t)=\frac{P_i(t-1){\epsilon }^T(t)}{1+\epsilon (t)P_i(t-1){\epsilon }^T(t)}$$

(40)

$$P_i(t)=[I-L_i(t)\epsilon (t)]P_i(t-1)$$

(41)

where $\widehat{\lambda }(t)$ and $\widehat{\eta }(t)$ represents the estimated parameters at time t. Notably, the parameters identified are $a_{21}$,$a_{22}$,$a_{23}$,$a_{24}$,$a_{29}$,$a_{210}$,$a_{31}$,$a_{61}$,$a_{1211}$,$a_{1213}$,$a_{1214}$,$a_{1311}$,$a_{1313}$,$a_{1314}$,$a_{61}$,$a_{1711}$,$a_{1714}$,$a_{1814}$,$n_2$ and $n_6$. Others are known.

3.3. Posting Process

Decomposing the combined parameters: the admittance parameters such as $J_{dr}$, $B_{dr}$, $K_{\mathrm{int}}$ can reflect the characteristics of the driver significantly. However, we only recognize the combined parameters until the last step. Therefore, the system of parameters of five variables are rewritten as follows:

$$\{\begin{array}{l}{a}_{{}_{21}}=-(\frac{K_{dr}-K_{st}+K_a-K_{\mathrm{e}}+1}{J_{dr}+J_{st}});a_{22}=-(\frac{B_{dr}+B_{st}-B_e}{J_{dr}+J_{st}})\hfill \\ a_{{}_{23}}=\frac{1}{J_{dr}+J_{st}};a_{31}=-\frac{W_c{B}_r}{I_r};a_{61}=-(W_c{K}_r-W_c{}^2{B}_r)\hfill \\ a_{1211}=\frac{I_{ped}-K_{tend}}{I_{seg}};a_{1213}=\frac{I_{ped}}{I_{seg}{B}_{con}};a_{1214}=\frac{K_{tend}}{I_{seg}{B}_{\mathrm{int}}}\hfill \\ a_{1313}=\frac{I_{ped}-K_{con}}{b_{con}};a_{1414}=-\frac{K_{\mathrm{int}}+K_{tend}}{b_{\mathrm{int}}}\hfill \\ a_{1711}=-K_f{K}_{tend};a_{1814}=\frac{K_{pos}-K_{vel}}{\tau b_{\mathrm{int}}}\hfill \end{array}$$

(42)

where $J_{dr}$, $B_{dr}$, $K_{dr}$, $B_r$, $K_r$, $K_{tend}$, $K_f$, $K_{vel}$, $K_{\mathrm{int}}$, $B_{\mathrm{int}}$, $K_{con}$, $B_{con}$ are unknown parameters, which can be calculated by other known vectors.

Furthermore, a few time-varying parameters of the driver’s NMS are difficult to be measured by physical sensors. The UKF observer based on the unscented transform (UT) theory and statistical linearization technique is implemented to observe the primary NMS parameters online, which reflect the evolution of driving skills [29].

4. Field Experiments

4.1. Experimental Facilities

The field experiments are conducted with Chongqing Jiaotong University’s “XinYuan” autonomous vehicle, as shown in Figure 4. All of the vehicle’s signals are collected by a V-BOX integrated device, a vehicle performance test system provided by the British Racelogic company. As shown in Figure 4, the steering wheel angle and torque is measured by a steering sensor mounted on the steering wheel. The pedal force and string position sensors are responsible for collecting the pedal force and stroke. The vehicle’s position and speed information are recorded by the Global Positioning System / Real-time kinematic (GPS/RTK). The vehicle’s state, such as $V_y$, ${\omega }_{\gamma }$, $e_y$, $\psi$, can be collected by the IMU inertial sensor.

In order to identify the parameters of the arm’s NMS characteristic, a 16-channel EMG measurement instrument supplied by the OT Bio company in Italy is utilized in the experiment. Its sampling rate, resolution rate and bandwidth are 1024 Hz, 16bit and 10 Hz–500 Hz, respectively. The Figure 5 illustrates that primary working muscle groups are monitored by EMG electrodes which are placed and oriented according to surface Electromyography for the non-invasive assessment of muscles (SENIAM) standards [30]. The working muscle groups affecting the driver’s steering ability are listed in

Table 3. Measured muscle groups.

Number	Muscle Name	Functionality
1	Sternocostal portion of Pectoralis major (PMA-S)	Steering right
2	Middle deltoid (DELT-M)	Steering right
3	Biceps brachii (BB)	Steering left
4	Lateral head of triceps brachii (TB-LA)	Steering right
5	Long head of triceps brachii (TB-L)	Steering left
6	Teres major (TM)	Steering left
7	Rectus femoris (RF)	Speed control
8	Patellar anterior (PA)	Speed control
9	Gastrocnemius muscle (GN)	Speed control
10	Soleus muscle (SL)	Speed control

[31,32].

4.2. Experimental Participants

With the agreement of our university’s academic committee and safety regulatory authorities, three subjects were recruited to take part in the experiments. All of them were male students in college, 26 years old with two years of driving experience. They were volunteers, taking part in the research with no financial reward. All the subjects had experiments performed three times, and the best one is chosen as the experimental data.

4.3. Experimental Scenario

The driving test field of Chongqing Jiaotong University is adopted as our experimental scenario, as shown in Figure 6. The red line represents the trajectory collected by GPS/RTK in the right lane. The complete road section consists of a 326 m double lane with four turning segments, four straight line segments and a ramp segment. The average driving speed on the experimental scenario is restricted to 30 km/h for safety.

4.4. Experimental Procedure

To enhance the confidence of the experimental results, each subject had 20 min of adaptive driving training. Furthermore, sampling frequency of EMG instruments and VBOX were adjusted to be as consistent as possible, both of which were 1024 Hz. Before the experiment, the maximum voluntary contraction (MVC) of each muscle working on steering wheel and pedal was calibrated in order to solve the muscle activation. During the experimental process, three subjects kept holding the steering wheel sensor at the “3–9 o’clock” position, and pressed the brake pedal sensor as shown in Figure 7. Each subject followed the trajectory three times, and the best instance was selected as the experimental data. Noteworthily, an emergency acceleration and braking scenario were designed in order to observe the fast response of the vehicle and muscle signals in S3 segments.

5. Model Verification with Experimental Results

5.1. Activation-to-Force Processing

Muscle force plays a vital role in the research of the driver’s NMS dynamic model. It cannot be measured directly, but can be obtained based on EMG signals. Therefore, a numbers activation-to-force model was proposed [26]. In this paper, the EMG signals are obtained by EMG analysis software “BioLab”. In order to facilitate the analysis of the result, the measurement noise is filtered through 20–400 Hz frequency, followed by average rectified value (RMS) processing. Finally, the linearized time-domain signal is obtained. Noteworthily, the following assumptions are presented: the driver’s muscle characteristic is linear; ignoring the influences of driving fatigue; the length of the tendon contraction is constant; the contribution of each work muscle can be quantified by different weights.

The muscle activation $a(t)$ is processed by a nonlinear function:

$$a(t)\approx u(t)=\frac{RMS(t)}{RM{S}_{\max }}$$

(43)

where $RMS(t)$ represents the RMS value of EMG at time t after normalization; $RM{S}_{\max }$ infers to the RMS value of EMG in isometric maximum voluntary contraction(MVC); $u(t)$ infers to the strength of muscle; $a(t)$ indicates the degree of muscle activation.

Ten muscle activations are calculated as shown in Figure 8. Evidently, the muscles of Biceps brachii (BB), lateral head of triceps brachii (TB-LA), pectoralis major (PMA-S) and long head of triceps brachii (TB-L) have higher activation which means more contributions for steering movement. Similarly, the rectus femoris (RF) muscle plays the most significant role during speed control.

5.2. Model Verification Results

Using the aforementioned identification algorithm with the experimental results, the key parameters of three subjects are identified in

Table 4. Identified parameter of NMS characteristics.

	Subject	Range	Initial	Sub1	Sub2	Sub3
Prameter
Jdr (kg/m2)		[0.03~0.14]	0.172	0.162	0.097	0.125
Bdr (Nm/rad)		[0.1~2.5]	1.032	0.989	0.91	0.932
Kdr (Nm/rad)		[20~90]	60.707	58.873	49.105	52.34
Br (Nm s/rad)		[0.5~2]	1.69	1.532	0.9	1.472
Kr (Nm/rad)		[2~30]	3.424	2.502	2.15	2.27
Ktend (Nm/rad)		[1799–3988]	2799	2771	2096	2459
Kf (Nm/rad)		[−0.81~1.78]	1.18	1.27	1.36	1.32
Kvel (Nm/rad)		[−0.40~48.6]	40.4	41.51	33.12	38.74
Kint (Nm/rad)		[191~505]	334	318	381	397
Bint (Nm/rad)		[16~22.2]	19.4	19.3	16.1	17.8
Kcon (Nm/rad)		[615~1143]	1033	1022	893	921
Bcon (Nm/rad)		[10.6~11.7]	11.2	11.32	10.87	10.91

. The original range of each parameter is derived from reference [33]. All NMS parameters are within a reasonable range, which indicates the model identification results are practicable. In particular, higher stiffness and damping coefficients indicate greater muscle force applied by the driver.

In order to prove the accuracy of the dynamic model and to observe the driver’s personalized behavior, the estimated value with the UKF observer of the first subject and actual value based on the experiment are illustrated in Figure 9, Figure 10 and Figure 11. As can be observed, the UKF estimated value has good consistency with the actual value.

5.3. Model Verification Analysis

In order to verify that the model can observe the driver’s NMS characteristics, the principle of covariance is utilized to analyze the correlation between the actual and observed value of Subject 1. From the statistical principle, the larger $\mathrm{var}(t)$ value represents the higher degree of data fusion between two groups.

$$\mathrm{var}(t)=\frac{{\mathsf{\Sigma }}_{t=1}^n(y_I(t)-{\overline{y}}_I)(y_s(t)-{\overline{y}}_s)}{n-1}$$

(44)

where, $y_I(t)$ and $y_s(t)$ represents the initial value and the estimated value which characterizing the subject at time t respectively; ${\overline{y}}_I$ and ${\overline{y}}_s$ refers to the mean value of the initial value and estimated value respectively.

As shown in

Table 5. Covariance analysis results between the actual value and Unscented Kalman Filter (UKF) estimated value.

No.	Parameter	Coefficient	No.	Parameter	Coefficient
1	Lateral speed	0.8621	9	Steering assistance torque	0.7156
2	Yaw rate	0.8246	10	Motion of leg	0.82323
3	Lateral offset	0.7653	11	Movement speed of leg	0.74421
4	Yaw	0.8956	12	Pedal torque	0.7021
5	Steering wheel angle	0.8663	13	Leg contraction torque	0.74421
6	Steering wheel angular rate	0.8421	14	Leg tendon motion	0.746
7	Arm active contraction torque	0.73421	15	Pedal position	0.716
8	Steering wheel torque	0.8646

, fifteen correlation coefficients between the actual value and the UKF estimated value are calculated with covariance analysis. From the analysis, all the coefficients greater than 0.7, which indicated that the parameter identification, as well as the UKF observer we proposed, have significant precision.

Furthermore, in order to prove the model’s practicability in reflecting the personalized driver characteristics, the actual value of Subject 1 is compared with three subjects’ estimated values. The covariance analysis results in Figure 12 is normalized for visualization. It can be observed that the histogram of Subject 1 has the highest covariance with the initial state vector. Finally, we can deduce that the HVSC dynamic model we proposed can represent the evolution of personal characteristics.

6. Conclusions

In this paper, a HVSC dynamic model composed by driver’s NMS model, 2DOF vehicle dynamic model and MPC controller is proposed for the shared control of AV. Specifically, the influence of the significant factors of the arms and legs, such as stiffness, damping and inertia are considered into the driver’s NMS model. Field experiments were conducted for the validation of the model. Based on the subsequent experimental results, the HLS algorithm and UKF observer are implemented to identify and estimate the important NMS parameters, respectively. Finally, the covariance analysis shows consistency with the identified parameter of the driver’s NMS characteristics. The results show that the model we proposed has a high level of accuracy with actual value. This paper will serve as a theoretical basis for the control strategy allocation between human and vehicle for L3 class AVs.