False cue influence on motion cue quality for 10 motion cueing algorithms

Abstract

Motion simulators are becoming increasingly popular for many applications in which human sensation is important to replicate and optimize target motions. For the emulation of the perceived human acceleration, motion cueing algorithms (MCAs) have been proposed in the literature that mimics the motion sensation by a combination of actual acceleration and tilted gravity effects, termed g-force or specific force. However, their relative performance has not yet been analyzed. This paper reviews existing families of MCAs and compares their performance for a simple offline S-shaped planar test trajectory featuring only lateral acceleration. The comparison is carried out both numerically using two previously published objective measures, the “performance indicator” of Pouliot, Gosselin, and Nahon, and the “good criterion” of Schmidt, as well as subjectively by a preliminary passenger rating on a real motion platform—Robocoaster testbed. The results show that (a) the novel optimizing MCA group exploits more effectively the workspace of the motion platform than the traditional MCA group for reducing false cue with small scale error and shape errors, (b) path-dependent tuning of MCA parameters may improve motion sensation, (c) average subjective ratings can be made to correlate well with the “good criterion” when expanded with a penalty for false angular velocity cues, and (d) the scale error of specific force seems to play the most important role to the evaluation of test subject on the motion cue quality. However, still a strong variance in subjective ratings was observed, making further research necessary.

Introduction

Driving simulation systems are useful and affordable system applying on both training and researching human driving behavior. The systems have first used the parallel robot named Stewart platforms to generate the supporting motion. For reducing these costs and, most notably, extending the motion workspace, new motion platforms based on serial robots have been recently developed, enlarging the applicability of the simulating process. One example of these serial-robot ride simulators is the KUKA Robocoaster, used in the entertainment industry as well as for research purposes. ¹ Common to all motion platforms, including serial robots, is that exact target motions cannot always be precisely reproduced due to the limited workspace. In these cases, the motion cueing algorithms (MCAs) play an important role in mincing the benefit motion sensations that can be reproduced by motion platforms in the restricted available workspace. MCAs divide the target acceleration into the high-and low-frequency parts. While the high-frequency (HF) parts are reproduced by the transient linear acceleration of the motion platform, the low-frequency (LF) parts of acceleration are simulated by the tilt coordination technique. The principle of the technique is to tilt the motion platform a suitable angle to treat the human vestibular system. However, the interaction between the threshold of human motion perception and the restricted available workspace of motion platforms often leads to false cues. The negative cues induce simulator sickness because the visual, vestibular, motion perception have no consistency information. For example, when perception thresholds for tilt-angle g-force emulations are exceeded, drivers can feel the rotational motions that do not exist in the target motion. Recently, various MCAs have still been developed for different motion platform, thus, for a specific simulation task, selecting the best MCA and tuning its parameters are still the open problems.

Existing MCAs employ different techniques to blend tilt-angle g-force and proper motion acceleration. At the origin of MCAs, classical MCA (or named washout filter) was proposed by Conrad and Schmidt ² consisting of a low-pass filter (LP) cross-coupled with a high-pass filter (HP), whose parameters were manually adjusted offline. Due to its simplification in tuning and deploying, this algorithm has been broadly used in commercial driving simulation systems. However, the ineffective exploitation of the available workspace is the weakness of the algorithm. Besides, the classical MCA generates false cues due to the linear characteristics of the washout filters. For better exploitation of the serial robot workspace, such as the KUKA Robocoaster, the classical MCA was modified into “cylindrical classical algorithm” by Giordano et al. ¹ The authors implemented the washout filters in a defined cylindrical coordinate by which the end-effector can move along both the linear and circular motion. Therefore, the cylindrical MCA can generate larger lateral accelerations than the classical MCA. In order to adjust parameters of washout filters in realtime, the “adaptive washout filter” was introduced.^3–5 The cost function of the adaptive method uses the steepest descent technique that is to minimize motion cue errors and to restrict linear movements of motion platform deviation of the filtered acceleration signal from the target signal. The adaptive MCAs decreases the false cues compared with the classical MCA, however, the stability of adaptive MCAs significantly depends on its parameters. ⁶ Further development has been to design MCAs that minimize a cost function involving the assumed perception error at the vestibular system and energy terms over the complete track. This family of methods, firstly introduced by Sivan et al., ⁷ solves the stochastic optimization problem to find higher-order washout filters for four couple motions. The Optimal MCAs were later improved by Nahon and Reid ⁶ and Telban et al. ⁸ who use different mathematical models of vestibular systems. Zywiol and Romano ⁹ considered on the tracking problem of simulated and target signal to propose a time backward optimal MCA. The algorithm solves the tracking problem with the linear model of tilt-coordination in offline mode to find the best exploitation of the available workspace for both translational and rotational motion. More recently, also novel model predictive control (MPC) was applied to the problem of motion cueing. The algorithm takes the various restriction into account, such as workspace limits, the motion perception thresholds (angular velocity, perception of linear acceleration), etc.^10–12 These algorithms show substantial improvements in motion sensation fidelity and are able to fill better the available workspace, but still, have some limitations in terms of computational efficiency and avoidance of false cues. Khusro et al. ¹³ proposed a nonlinear MPC solution that uses the Stewart platform’s nonlinear kinematics into the MPC algorithm to improve cueing fidelity and maximize workspace. The author also employed the adaptive weights tuning to smooth the platform’s movement toward its physical limits. ¹³ Katliar et al.^14,15 developed an MPC-based MCA for a Cable-Robot-based motion simulator that includes motion platform actuation. Their key conclusion was that an MPC-based MCA with a complex model can be executed in real time with the right software and numerical approaches. Katliar et al. ¹⁶ used a model predictive controller (MPC) for real-time control of a motion simulator based on an 8-degree-of-freedom serial robot. By moving the human participant sitting within the cabin positioned at the end effector, the controller aims to properly duplicate six reference signals (accelerations and angular velocities in the body frame of reference) acquired from a simulated or actual vehicle. Besides, the use of a nonuniform weighting strategy to stabilize the MCA using MPC with a short prediction horizon and improved weighting modification is proposed in the study of Mohammadi et al. ¹⁷ The author mentioned that reduced prediction horizons are beneficial because they reduce computing strain, but they also move the system toward instability. Moreover, Asadu et al. ¹⁸ proposed a method based on compensators and generic algorithm that has the goal of minimizing the fitness function, which can lead to greater fidelity motions and better solutions. The following factors are considered when creating the fitness function such as minimizing the translational and rotational human sensation error between the real and simulator drivers; maximizing the correlation coefficient; minimizing the angle and linear displacements; minimizing the acceleration, velocity, and displacement of the platform; and minimizing sensation error fluctuations.

The objective of the present work is (1) to review currently available MCAs for the families of classical, adaptive, optimal control, and model-predictive control algorithms in offline mode; (2) to numerically compare their performance for a simple planar S-shaped test ride using two existing objective measures for perception evaluation, the “performance indicator” of Pouliot et al., ¹⁹ and the “good criterion” of Fischer ²⁰ ; and (3) to correlate the thus obtained objective measures with preliminary subjective ratings on a real simulator ride.

The rest of the paper is structured as follows. First, the testbed is described in Section 2. Section 3 then reviews briefly all existing MCAs, whose perception quality on the test ride is evaluated numerically in Section 4 and by preliminary subjective ratings in Section 5. The results finally are discussed in the conclusions section.

Driving simulator based on KUKA robocoaster KR500

The motion platform based on KUKA robocoaster is applied for a driving simulator due to its flexibility in rotational motions that are restricted in the steward motion platform. For example, the motion platform in the Laboratory for Mechanics and Robotics (LMR) at University of Duisburg-Essen uses a KUKA KR500/1 TÜV Robocoaster robot, shown in Figure 1. The motion platform can lift driver seating on a Maurer Söhne rollercoaster seat at the robot end effector. Besides, the head-mounted display unit—Visette45—provides the virtual reality of a 3D environment to the driver. The trajectory of the motion platform is generated by motion cueing algorithms and then transferred to the robot’ controller. The virtual reality environment is simultaneously computed to the trajectory of the robot in the computer visualizing environment, then is delivered in stereo 3D to the head-mounted display. The physical limits of the robot’s joints are shown in Table 1. When compared to the Stewart platform, the KUKA robocoaster can reproduce a larger lateral motion by utilizing its enormous circular motion due to the large rotating space of the base joint (joint 1). On the other hand, the simulation of vertical motion is still restricted in a range of 0.127 m.

Figure 1.

Architecture of the LMR driving simulator (base on Pham et al. ²¹ ).

Table 1.

Specification of the KUKA Robotcoaster.

Joint	Maximum angular acceleration (rad/s2)	Maximum angularvelocity (°/s)	Angular hard-limit
			Min (°)	Max (°)
1	1.7	69	−167	167
2	1.2	69	−96	−15
3	2.2	69	−28	17
4	0.7	77	−350	350
5	1.6	76	−57.5	57.5
6	1.1	120	−350	350

Review of existing MCAs

There exists a great variety of motion cueing algorithms (MCAs) in the current literature. They are mainly devoted to mimic motion perception by a superposition of actual passenger acceleration with tilted gravity vector effects, termed the specific force. In a driving simulator, the motion cues are dominantly sensed by the vestibular system that perceives translational and rotational motion represented by specific force and angular velocity, respectively. ²² This section reviews the principal MCA families with particular interest on the transversal accelerations ${{}^P\mathit{a}}_v=[0,a_{\mathrm{vy}},0{]}^T$ in body-fixed coordinates, where $x,y,z$ are the “front,”“lateral,” and “up” coordinates, respectively, and the target angular velocity is ${{}^P\mathit{\omega }}_v=[{\omega }_{\mathrm{vx}},{\omega }_{\mathrm{vy}},{\omega }_{\mathrm{vz}}{]}^T=[0,0,0{]}^T$ . The rotation about the front axis is measured by the roll angle ${\phi }_v$ , and the linear acceleration is measured at a reference point “V” in the cabin, which is assumed to be located at the passenger’s head “P.”

Classical MCA with cylindrical coordinates (CL)

This algorithm (Figure 2), introduced by Giordano et al., ¹ applies the standard classical washout algorithm to cylindrical coordinates. Hereby, absolute Cartesian coordinates $\mathit{p}=[X,Y,Z{]}^T$ of the cabin (Z is “up”) are transformed to cylindrical coordinates according to equation (1).

$$\mathit{\xi }=[R,\alpha ,Z];\mathrm{with}R=\sqrt{X^2+Y^2};\tan (\alpha )=\frac{Y}{X}.$$

(1)

Figure 2.

Control diagram of classical MCA in cylindrical coordinate. ¹

At the robot base, a vertically rotating frame $K_W$ is attached such that the cabin reference point always lies in its $x/z$ plane. The cabin orientation is measured with respect to this frame by pitch and roll rotations ${\theta }_W,{\phi }_W$ , in that order, collected in the vector ${\mathit{\beta }}_W$ . The target acceleration ${{}^P\mathit{a}}_v$ is first transformed to $K_W$ by a rotation matrix ${{}^{W}R}_P$ and then to cylindrical accelerations $\stackrel{··}{\mathbf{\xi }}$ by the linear mapping $C(\xi )=\mathrm{diag}\{1,\frac{1}{R},1\}$ . These are then fed into the three channels of linear high-pass filters (equation (2)) yielding high-frequency components ${\stackrel{··}{\mathbf{\xi }}}_{\mathrm{HP}}$ . Due to the circular motion, centripetal and Coriolis accelerations $-R{\stackrel{·}{\alpha }}_{\mathrm{HP}}^2$ and $2\stackrel{·}{R}{\stackrel{·}{\alpha }}_{\mathrm{HP}}$ , respectively, are generated which must be compensated by cabin pitch and roll tilt angle contributions ${\theta }_T=-\arcsin \frac{{{}^{P}a}_{\mathrm{Tx}}}{g}$ , ${\phi }_T=\arcsin \frac{{{}^{P}a}_{\mathrm{Ty}}}{g\cos {\theta }_T}$ , respectively, collected in vector ${\beta }_{\mathrm{IN}}$ . These are added to the result of the low-pass filter $\mathrm{LP}(s)=1-\mathrm{HPa}(s)$ applied to the target acceleration. The final cabin pitch/roll angles ${\mathit{\beta }}_W$ are assumed to be the sum of ${\mathit{\beta }}_T$ and the pitch/roll angles ${\mathit{\beta }}_{\mathrm{HP}}$ from the time integral of the high-pass filter (equation (3)) applied to the target angular velocity transformed to cabin coordinates. Prepending the rotation of the frame $K_W$ to the pitch/roll angles ${\beta }_W$ as a yaw angle gives the global orientation ${\beta }_s$ of the cabin as yaw/pitch/roll angles (in that order). The filter parameters for all channels used in our comparison are shown in Table 2, where $K$ is the proportional gain, $\zeta$ is the damping ratio, ${\omega }_n$ is the natural frequency and ${\omega }_b$ is the break frequency.

$$\mathrm{HPa}\left(s\right)=\frac{K{s}^2}{s^2+2\zeta {\omega }_{n}s+{\omega }_n^2}·\frac{s}{s+{\omega }_b},$$

(2)

$$\mathrm{HP}\omega (s)=K{s}^2/(s^2+2\zeta {\omega }_{n}s+{\omega }_n^2)$$

(3)

Table 2.

Classical algorithm’s parameters. ¹

	$K$	$\zeta$	ωn (s−1)	ωb (s−1)
$\{R;\alpha \}$	$\{1;1\}$	$\{0.7;1\}$	$\{0.5;0.5\}$	$\{1,2;1.2\}$
$\{z;\omega \}$	$\{1;1\}$	$\{0.7;0.7\}$	$\{0.6;1\}$	$\{20;-\}$

Adaptive MCA with cylindrical coordinates (AD)

The algorithm (Figure 3) is developed from the classical MCA with cylindrical coordinates by using a steepest descent method in an adaptive block “ $\mathrm{ADAPT}$ ” to adapt the parameter set $\{K$ , ${\omega }_n$ , $\zeta \}$ of the high-pass filter $H{P}_2=\frac{K{s}^2}{s^2+2\zeta {\omega }_{n}s+{\omega }_n^2}$ in real time to minimizing the following cost function (equation (4)).

$$\mathit{J}=\{\frac{1}{2}[w_a({\mathit{\xi }}_F-{\mathit{\xi }}_{\mathrm{HP}}{)}^2+w_v{\mathit{\xi }}_{\mathrm{HP}}^2+w_s{\mathit{\xi }}_{\mathrm{HP}}^2+w_k(K-K_0{)}^2+w_{\zeta }(\zeta -{\zeta }_0{)}^2+w_{\omega }({\omega }_n-{\omega }_{n0}{)}^2]\}$$

(4)

Here, ${\mathit{\xi }}_F$ is the filtered acceleration through $H{P}_1=\frac{s}{s+{\omega }_b}$ , $K_0,{\zeta }_0,{\omega }_{n0}$ are the reference anchor values of the filter parameters, and $w_{(\cdots )}$ are weighting parameters (Table 3), which penalize the deviation between the target and the simulated acceleration, and restrain the cabin’s movement.

Figure 3.

Control diagram of adaptive MCA in cylindrical coordinates. ⁵

Table 3.

Parameters of adaptive MCA in cylindrical coordinates. ⁵

Initial parameters			Weighting parameters
$K_0$	ωn 0 (s−1)	${\zeta }_0$	ωb (s−1)	$w_k$	$w_a$	$w_{\zeta }$	$w_{\omega }$	$w_v$	$w_s$
1	0.5	0.7	1.2	1	1	1	0.01	1	1

Optimal washout MCAs (OpS, OpNR, Tel-YoMe, and Tel)

The problem was first raised by Sivan et al. ⁷ The authors defined the linear quadratic optimal control for tilt coordination problem to find higher-order filters $W(s)$ (Figure 4). The filters minimize the cost function (equation (5)), and the state space model (equation (6)).

$$\mathit{J}=E\left\{{\int }_{t_0}^{t_f}({\mathit{e}}^{\mathit{T}}\mathit{Qe}+{\mathit{u}}^{\mathit{T}}\mathit{Ru}+{{\mathit{x}}_{\mathit{c}}}^{\mathit{T}}{\mathit{R}}_{\mathit{c}}{\mathit{x}}_{\mathit{c}})\mathrm{dt}\right\}\stackrel{!}{=}\min$$

(5)

$$\{\begin{array}{c}\stackrel{·}{\mathit{x}}=\mathit{Ax}+\mathit{Bu}+\mathit{Hn}\hfill \\ \multicolumn{1}{c}{\mathit{y}=\mathit{Cx}+\mathit{Du}}\end{array}$$

(6)

Here, $E\{\}$ is the mean estimator, and $[t_0,t_f]$ is the time-simulation time. The output vector $\mathit{y}={[\mathit{e},{\mathit{x}}_c]}^T$ contains the sensation error $\mathit{e}(t)=[\mathrm{\Delta }{\hat{f}}_y,\mathrm{\Delta }{\hat{\omega }}_x{]}^T$ of the specific force and angular velocity, respectively, and the ${\mathit{x}}_c$ simulator state variables. It is important to note that the state variable $x$ , system inputs $u$ and $x_c$ , as well as the weighting matrices $\mathit{R}=\mathrm{diag}\{r_1,r_2\}$ , $\mathit{Q}=\mathrm{diag}\{q_1,q_2\}$ , and ${\mathit{R}}_c=\mathrm{diag}\{r_{c1},\dots ,r_{c4}\}$ , differ between methods and are given in Tables 4 and 5, respectively. For example, in the basic approach Sivan (equation (6)) defined the state variable: $\mathit{x}={[x_v,x_s,x_n]}^T$ , where the state variables in the vehicle ${\mathit{x}}_v={\left[{\mathit{x}}_v^{oto},{\mathit{x}}_v^{ssc},\int {\int }^{}{a}_{vy},{\int }^{}{a}_{vy}\right]}^T$ with otolith state variable ${\mathit{x}}_v^{\mathrm{oto}}=[{\hat{f}}_{\mathrm{vy}}-G_O{B}_o{f}_{\mathrm{vy}}]$ , and semi-circular state variable ${\mathit{x}}_v^{ssc}={\left[{\widehat{\omega }}_{vx}-G_S{b}_S{\phi }_v,{\int }^{}{\widehat{\omega }}_{vx}\right]}^T$ . Similarly, in term of the simulator state variables ${\mathit{x}}_s$ is defined as ${\mathit{x}}_v$ (with subscript “v” is replaced by “s”). Note that ${\mathit{x}}_n=[a_{\mathrm{vy}},{\phi }_v]$ and ${\mathit{H}}_n$ are the state of white noise filter and state-space noise, respectively. The matrices in the first equation of the dynamic model in (equation (6)) are represented as equation (7).

$$\mathit{A}=\left[\begin{array}{ccc}\overline{\mathit{A}}\hfill & 0\hfill & {\overline{\mathit{B}}}_{\mathit{v}}\hfill \\ \multicolumn{1}{c}{0}& \overline{\mathit{A}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & {\mathit{A}}_{\mathit{n}}\hfill \end{array}\right],\mathit{B}=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{\overline{\mathit{B}}}\\ \multicolumn{1}{c}{0}\end{array}\right],\mathbf{and}\mathit{H}=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{{\mathit{H}}_{\mathit{n}}}\end{array}\right]$$

(7)

Figure 4.

Control diagram of optimal washout algorithms. ⁸

Table 4.

Dynamic states applied in the model of optimal problems.

	${\mathit{x}}_c$				$\mathit{u}$		$\mathit{\gamma }$
	${\int \int \int a_{sy}d{t}^3}^{\$}$	$\int {\int }^{}{a}_{sy}d{t}^2$	${\int }^{}{a}_{sy}dt$	${\phi }_s$	$u_1$	$u_2$	γ 1 (s−1)	γ 2 (s−1)
OpS 7	−	+	+	−	$a_{\mathrm{sy}}$	${\phi }_s$	1	0.1
OpRN 4	+	+	+	−	${\phi }_s$	$a_{\mathrm{sy}}$	0.2	1
Tel-YoMe 8	+	+	+	+	${\stackrel{·}{\phi }}_s$	$a_{\mathrm{sy}}$	1	4* $\pi$
Tel 8	+	+	+	+	${\stackrel{·}{\phi }}_s$	$a_{\mathrm{sy}}$	1	$\pi$

“−”: the abandoned parameters; “+”: the used parameters.

Table 5.

Weighting parameters of the optimal washout algorithms.

	$q_1$	$q_2$	$r_1$	$r_2$	$r{c}_1$	$r{c}_2$	$r{c}_3$	$r{c}_4$
OpS 7	0.707	0.707	–	0.999	–	$r{c}_S$	$r{c}_S$	–
OpRN 4	0.447	0.894	0.9	$r{c}_R$	$r{c}_R$	$r{c}_R$	$r{c}_R$	–
Tel-YoMe 8	1	10	1	1	8	4	1	250
Tel 8	1	10	1	1	8	4	1	250

$r{c}_S=0.00577$ and $r{c}_R=0.158$ .

while, the second equation has the matrices as equation (8).

$$\mathit{C}=\left[\begin{array}{ccc}-\overline{\mathit{C}}\hfill & \overline{\mathit{C}}\hfill & -\overline{\mathit{D}}\hfill \\ \multicolumn{1}{c}{0}& {\mathit{C}}_{\mathit{c}}\hfill & 0\hfill \end{array}\right],\mathit{D}=\left[\begin{array}{c}\overline{\mathit{D}}\hfill \\ \multicolumn{1}{c}{0}\end{array}\right].$$

(8)

Moreover, the state matrix ${\mathit{A}}_{\mathit{n}}=\mathrm{diag}\{-{\gamma }_1,-{\gamma }_2\}$ and input matrix ${\mathit{B}}_{\mathit{n}}=\mathrm{diag}\{{\gamma }_1,{\gamma }_2\}$ are used for modeling white noise transmission, ${\mathit{C}}_{\mathit{c}}=[0,I]$ , while the matrices $\overline{\mathit{A}}$ , $\overline{\mathit{B}}$ , $\overline{\mathit{C}}$ , $\overline{\mathit{D}}$ integrate the sub-matrices of vestibular dynamic model as $\tilde{\mathit{A}}$ , $\tilde{\mathit{B}}$ , $\tilde{\mathit{C}}$ , $\tilde{\mathit{D}}$ depending on the method. For example, the basic approach of Sivan et al. ⁷ (OpS) uses equation (9).

$$\overline{\mathit{A}}=\left[\begin{array}{cc}\tilde{\mathit{A}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\mathit{A}}_{\mathit{c}}\hfill \end{array}\right],\overline{\mathit{B}}=\left[\begin{array}{c}\tilde{\mathit{B}}\hfill \\ \multicolumn{1}{c}{{\mathit{B}}_{\mathit{c}}}\end{array}\right],\overline{\mathit{C}}=\left[\begin{array}{c}\tilde{\mathit{C}}\hfill \\ \multicolumn{1}{c}{0}\end{array}\right],\overline{\mathit{D}}=\tilde{\mathit{D}},$$

(9)

with equation (10).

$$\begin{array}{c}\tilde{\mathit{A}}=\left[\begin{array}{ccc}-B_0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -a_S\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& -b_S\hfill & 0\hfill \end{array}\right],\tilde{\mathit{B}}=\left[\begin{array}{cc}{G}_O(A_0-B_0)\hfill & -G_{O}g(A_0-B_0)\hfill \\ \multicolumn{1}{c}{0}& -a_S{b}_S{G}_S\hfill \\ \multicolumn{1}{c}{0}& -b_S^2{G}_S\hfill \end{array}\right],\hfill \\ \tilde{\mathit{C}}=\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \end{array}\right],\tilde{\mathit{D}}=\left[\begin{array}{cc}{G}_O\hfill & -G_{O}g\hfill \\ \multicolumn{1}{c}{0}& G_S{b}_S\hfill \end{array}\right],{\mathit{A}}_{\mathit{c}}=\left[\begin{array}{cc}0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right],{\mathit{B}}_{\mathit{c}}=\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill \end{array}\right].\hfill \end{array}$$

(10)

Here, $A_0,B_0,G_O,G_S$ are the physical parameters of the dynamic models of vestibular system (Tables 6 and 7), and $a_S=({\tau }_1+{\tau }_2)/{\tau }_1{\tau }_2$ and $b_S=1/{\tau }_1{\tau }_2$ .

Table 6.

The otolith parameters (sway direction).

	A 0 (s−1)	B 0 (s−1)	B 1 (s−1)	G 0 (s2/m)	δ 0 (m/s2)
OpS 7	0.076	0.19	−	2.16	0.47
OpRN 4	0.076	0.19	−	5.86	0.17
Tel-YoMe 8	0.076	0.19	1.5	8.82	0.17
Tel 8	0.1	0.2	62.5	294.12	0.17

“−”: the abandoned parameters.

Table 7.

Semicircular parameters (roll angle).

	τ 1 (s)	τ 1 (s) ${\tau }_1(s)$	τa (s)	τl (s)	Gs (s2/rad)	δ s (°/s)
OpS 7	5.9	0.003	–	–	233	$1.45$
OpRN 4	6.1	0.1	–	–	118.55	3
Tel-YoMe 8	5.73	0.005	80	0.06	28.6479	2
Tel 8	5.73	0.005	80	0.06	28.6479	2

Telban et al. ⁸ (Tel) extended these models to higher-order filters

$$\mathrm{scc}(s)=G_S\frac{{\tau }_a{\tau }_1{s}^2(1+{\tau }_{l}s)}{({\tau }_{1}s+1)({\tau }_{2}s+1)({\tau }_{a}s+1)}$$

(11)

and

$$\mathrm{oto}\left(s\right)=G_O\frac{(s+A_0)}{(s+B_0)(s+B_1)},$$

(12)

for the semicircular and otolith system, also using a Young-Meiry vestibular model (Tel-YoMe). An intermediate filter complexity was introduced by Nahon and Reid ⁶ (OpNR). Tables 6 and 7 summarize the parameters used in these optimal washout signals, where missing parameters mean that this order is not present in the transfer function. The scalar elements $W_{\mathrm{ij}}$ of the washout filter in Figure 4 are then determined in offline mode as the components of the $2\times 2$ matrix

$$\mathit{W}(s)=-{(\mathit{I}+{\mathit{F}}_2{(\mathit{sI}-{\mathit{A}}_{\mathit{s}})}^{-1}{\mathit{B}}_{\mathit{s}})}^{-1}·({\mathit{F}}_1{(\mathit{sI}-{\mathit{A}}_{\mathit{v}})}^{-1}{\mathit{B}}_{\mathit{v}}+{\mathit{F}}_3),$$

(13)

where ${\mathit{F}}_{\mathit{i}},i=1,\dots ,3$ are the submatrices corresponding to the subvectors ${\mathit{x}}_{\mathit{v}},{\mathit{x}}_{\mathit{s}},{\mathit{x}}_{\mathit{n}}$ of the feedback matrix in the expression $\mathit{u}=-\mathit{Fx}$ obtained as

$$\mathit{F}=[{\mathit{F}}_1,{\mathit{F}}_2,{\mathit{F}}_3]={\mathit{R}}_2^{-1}({\mathit{B}}^{\mathit{T}}\mathit{P}+{\mathit{R}}_{12}^{\mathit{T}}),{\mathit{R}}_1={\mathit{C}}^{\mathit{T}}\mathit{QC},{\mathit{R}}_{12}={\mathit{C}}^{\mathit{T}}\mathit{QD},{\mathit{R}}_2=\mathit{R}+{\mathit{D}}^{\mathit{T}}\mathit{QD}.$$

(14)

Note that the solution of the algebraic Riccati equation P is the unique and non-negative.

$$\mathit{P}(\mathit{A}-{\mathit{BR}}_2^{-1}{\mathit{R}}_{12}^{\mathit{T}})+{(\mathit{A}-{\mathit{BR}}_2^{-1}{\mathit{R}}_{12}^{\mathit{T}})}^{\mathit{T}}\mathit{P}+{\mathit{R}}_1-{\mathit{R}}_{12}{\mathit{R}}_2^{-1}{\mathit{R}}_{12}^{\mathit{T}}-{\mathit{PBR}}_2^{-1}{\mathit{B}}^{\mathit{T}}\mathit{P}=0$$

(15)

The washout filter works by first transforming the target acceleration $a_{\mathrm{vy}}$ through the block “ $\mathrm{TF}$ ” to inertial coordinates and then filtering it through $W_{11}$ to yield the simulator acceleration $a_{\mathrm{sy}}$ that passes through the integrator block to produce the simulator movement $S_y$ . Note that, the target acceleration $a_{\mathrm{vy}}$ is first filtered by $W_{12}$ —tilted coordination filter—and then the filtered part is summed with the filtered roll target angle through $W_{22}$ to produce the roll tilt angle ${\phi }_s$ . The cross-filter $W_{21}$ commonly has an insignificant gain, thus it can be neglected.

Offline optimal control algorithms (ZyRo and ZyRo*)

The algorithm was first introduced by Zywiol and Romano, ⁹ solves the tilt coordination as an optimal output tracking problem, whose cost functional (equation (16)) is minimized by the appropriate time-variant input vector (controlling vector) $\mathit{u}=[u_1,u_2{]}^T$ with $u_1$ representing for specific force at the pilot’s head, and $u_2$ representing for angular velocity. The input components $u_1$ , $u_2$ are passed through low-pass filter (PT1) to avoid the promptly unpractical change of the specific force and the angular velocity, respectively. The objective is that the output signal $\mathit{y}=[f_{\mathrm{sy}},S_y{]}^T$ containing the specific force and the lateral displacement of the cabin tracks the target output vector ${\mathit{y}}_{\mathit{ref}}=[f_{\mathrm{vy}},0{]}^T$ as close as possible while “pulling” the cabin to the home configuration. Let $\mathit{x}=[f_{\mathrm{sy}},{\stackrel{·}{\phi }}_s,{\phi }_s,{\stackrel{·}{S}}_y,S_y]T$ be the state variables and $\mathit{A},\mathit{B},\mathit{C}$ be the corresponding system matrices of the control system. The control input vector can be computed as equation (17).

$$\mathit{J}={\int }_{t_0}^{t_f}\left\{{\left[\mathit{y}-{\mathit{y}}_{\mathit{r}\mathit{e}\mathit{f}}\right]}^{\mathit{T}}\mathit{Q}\left[\mathit{y}-{\mathit{y}}_{\mathit{r}\mathit{e}\mathit{f}}\right]+{\mathit{u}}^T\mathit{R}\mathit{u}\right\}dt$$

(16)

$$\mathit{u}=-{\mathit{R}}^{-1}\left[{\mathit{B}}^{\mathit{T}}\mathit{U}\mathit{x}+{\mathit{B}}^{\mathit{T}}\overline{\mathit{S}}\right],$$

(17)

where matrix U, which is the positive semi-definite, is the solution of the following differential Riccati equation (equation (18)).

$$\stackrel{·}{\mathit{U}}=-\mathit{UA}-{\mathit{A}}^T\mathit{U}-{\mathit{C}}^T\mathit{QC}+\mathit{UB}{\mathit{R}}^{-1}{\mathit{B}}^T\mathit{U}$$

(18)

Note that $\overline{S}$ in (equation (17)) is the solution of co-state equations, and it is solved in backward time as in the equation (19).

$$\stackrel{.}{\overline{\mathit{S}}}=-[{\mathit{A}}^{\mathit{T}}-\mathit{UB}{\mathit{R}}^{-1}{\mathit{B}}^{\mathit{T}}]\overline{\mathit{S}}+{\mathit{C}}^{\mathit{T}}\mathit{Q}{\mathit{y}}_{\mathit{ref}}$$

(19)

The components of weighting matrices $\mathit{Q}=\mathrm{diag}\{q_1,q_2\}$ , $\mathit{R}=\mathrm{diag}\{r_1,r_2\}$ are tuned for tracking purposes of specific force and returning to the home position. On the other hand, the first break frequencies of low-pass filters (PT1) $c$ and $\gamma$ from Figure 5 are tuned for smooth response of angular velocity and specific force. The parameters are summarized in Table 8 for the original paper of Zywiol and Romano ⁹ (termed “ZyRo”) as well as for an alternative selection of values tuned for the motion platform based on KUKA robocoaster.

Figure 5.

Control diagrams of Zywiol Romano optimal algorithm. ⁹

Table 8.

Parameters of ZyRo algorithm.

	γ (s−1)	c (s−1)	$q_1$	$q_2$	$r_1$	$r_2$
ZyRo [?]	20	10	1	0.005	4	1
ZyRo*	20	100	199	0.4	0.01	199

Model predictive control algorithms (MPC and MPC*)

MPC algorithms base on model predictive control for the tilt coordination problem. The MPC is the multivariable optimization method that generates the best future behavior of a plant output by observing possible responses in a future fixed interval (prediction horizon). The possible responses are computed with corresponding parametrized future inputs over the prediction horizon while respecting the defined constraints. Therefore, the MPC algorithms can deal with both simulator limitations and human motion perception thresholds. Moreover, the algorithm can improve the motion cueing quality. There exist implicit MPC ¹⁰ based on an online optimization technique with Quadratic Programing (QP) as well as explicit MPCs ¹² based on an offline optimization technique.

Implicit MPC motion cueing algorithms

With the aim of reproducing the target specific force while keeping the simulator cabin inside the simulator workspace, Augusto and Loureiro ¹⁰ designed an implicit MPC algorithm by the minimization of the cost functional (equation (20)).

$$\mathit{J}(\mathrm{\Delta }\mathit{u})={[\mathit{y}-{\mathit{y}}_{\mathit{ref}}]}^{\mathit{T}}\mathit{Q}[\mathit{y}-{\mathit{y}}_{\mathit{ref}}]+{\mathit{u}}^{\mathit{T}}\mathit{Su}+\mathrm{\Delta }{\mathit{u}}^{\mathit{T}}\mathit{R}\mathrm{\Delta }\mathit{u},$$

(20)

where $\mathit{y}$ , ${\mathit{y}}_{\mathit{ref}}$ , $\mathit{u}$ , and $\mathbf{\Delta }\mathit{u}$ are the output vector, its reference values, the input vector, and its increment, respectively. Besides, $\mathit{Q}=\mathrm{diag}\{q_1,\dots ,q_n\}$ , $\mathit{S}=\mathrm{diag}\{s_1,\dots ,s_n\}$ , and $\mathit{R}=\mathrm{diag}\{r_1,\dots ,r_n\}$ are weighting matrices. Augusto and Loureiro ¹⁰ developed the MPC algorithm for four operating models: pitch/surge, roll/sway, heave, and yaw. Here, the lateral mode (roll/sway), which will be termed MPC in the following (Figure 6), are considered. The control vector consists of three inputs $\mathit{u}=[{\stackrel{·}{\phi }}_{\mathrm{as}},{\stackrel{·}{\phi }}_r,a_{\mathrm{sy}}{]}^T$ denoting the tilt roll rate for simulating specific force, the actual vehicle roll rate, and the simulated acceleration, respectively. For the output tracking problem, the solution of the MPC is the optimized control signal $u$ such that the output vector $\mathit{y}=[{\hat{\omega }}_{\mathrm{sx}},{\phi }_s,{\hat{f}}_{\mathrm{sy}},S_y,v_y{]}^T$ comprising of the sensation of roll angular velocity, the cabin roll angle, the lateral sensation of specific force, the velocity, and position of the cabin, respectively, tracks the corresponding reference signals ${\mathit{y}}_{\mathit{ref}}=[{\hat{\omega }}_{\mathrm{vx}},0,{\hat{f}}_{\mathrm{sy}},0,0{]}^T$ . The matrix $\mathit{H}$ in Figure 6 maps the tilt angle ${\phi }_s$ and simulated acceleration $a_{\mathrm{sy}}$ to the simulated specific force, $f_{\mathrm{sy}}$ , $\mathit{H}=[-g,1]$ . While applying the MPC algorithm with original parameters (Table 9 with $p,m$ are the number of time step for horizontal and control prediction, respectively), it turned out that the motion had a tendency to dephase with respect to the target specific force and to undulate in the angular velocity. Note that an unstable, jerky motion of the cabin could exist if the weighting matrices are selected to achieve the desired specific force fidelity. To avoid these problems, an alternative MPC algorithm termed MPC* (Figure 7) are proposed. The MPC has two inputs and insert, as in the optimal control algorithm, ⁹ two first-order low-pass filters in both the translational and rotational input channel (Figure 7, Table 9) to smooth the input signals. The break frequency of two low-pass filters were selected as in ZyRo* algorithm. Likewise, the output vector was reduced to $\mathit{y}=[{\hat{\omega }}_{\mathrm{sy}},S_y,{\hat{f}}_{\mathrm{sy}}{]}^T$ with reference values ${\mathit{y}}_{\mathit{ref}}=[{\hat{\omega }}_{\mathrm{vy}},0,{\hat{f}}_{\mathrm{vy}}{]}^T$ . Moreover, while the original MPC approach uses the extended linearized vestibular model proposed in Reid and Nahon, ⁴ we used the reduced model proposed therein without with additional time parameters $B_1=1.5(1/s)$ and ${\tau }_a=30(s)$ (Tables 6 and 7) to reduce computational effort. Besides, the horizon $p$ and control prediction $m$ parameters were chosen larger than those of original MPC to improve the tracking quality.

Figure 6.

Original control system of MPC algorithm.

Table 9.

Weighting values of MPC.

	$s_i$	$r_i$	$q_i$	$p$	$m$
MPC 4	$1|20|1$	$1|8|0.01$	$10000|20|10|10|15$	18	6
MPC*	$30|0.1$	$9000|0.1$	$50|5|99$	400	200

Figure 7.

Control system of MPC* algorithm.

Explicit MPC motion cueing algorithm

Based on the multi-parameter programing toolbox (MPT), developed by the Automatic Control Laboratory of ETH, Zürich, and the solution for the tracking problem, Fang and Kemeny designed the MCA based on an explicit MPC for real-time driving simulators (Figure 8), in which the optimized input vector minimizes the cost function (equation (21)) in the prediction horizon $N$ . Here, ${\mathit{u}}_{\mathit{i}}$ , ${\mathit{x}}_{\mathit{i}}$ , and ${\mathit{y}}_{\mathit{i}}$ denote the ith predicted input vector, predicted state vector, and predicted output vector, respectively; ${\mathit{x}}_{\mathit{N}}$ represents the final state vector and ${\mathit{y}}_{\mathit{ref}}{,}_{\mathit{i}}$ stand for the reference vector; the weighting matrices are diagonal matrices such as ${\mathit{Q}}_{\mathit{N}}=\mathrm{diag}\{q_1,\dots ,q_n\}$ , ${\mathit{Q}}_{\mathit{x}}=\mathrm{diag}\{q_1^x,\dots ,q_n^x\}$ , $\mathit{R}=\mathrm{diag}\{r_1,\dots ,r_n\}$ , and ${\mathit{Q}}_{\mathit{y}}=\mathrm{diag}\{q_1^y,\dots ,q_n^y\}$ .

$$\begin{array}{c}{\mathit{J}}_{\mathit{N}}({\mathit{x}}_{\mathit{k}})=\underset{{\mathit{u}}_0,{\mathit{u}}_1,...,{\mathit{u}}_{\mathit{N}-1}}{\mathbf{min}}{{\mathit{x}}_{\mathit{N}}}^{\mathit{T}}{\mathit{Q}}_{\mathit{N}}{\mathit{x}}_{\mathit{N}}+\sum _{\mathit{i}=0}^{\mathit{N}-1}{\mathit{u}}_{\mathit{i}}^{\mathit{T}}\mathit{R}{\mathit{u}}_{\mathit{i}}\hfill \\ +\sum _{\mathit{i}=0}^{\mathit{N}-1}{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{Q}{\mathit{x}}_{\mathit{i}}+\sum _{\mathit{i}=0}^{\mathit{N}-1}{({\mathit{y}}_{\mathit{i}}-{\mathit{y}}_{\mathit{ref}}{,}_{\mathit{i}})}^{\mathit{T}}{\mathit{Q}}_{\mathit{y}}({\mathit{y}}_{\mathit{i}}-{\mathit{y}}_{\mathit{ref}}{,}_{\mathit{i}})\hfill \end{array}$$

(21)

For smoothing motion in the available simulator’s workspace, Fang and Kemeny used a braking control law which avoids the jerky stop at the limit position (equation (24)). Starting from state-space model (equations (22) and (23))

$${\mathit{x}}_{\mathit{k}}=\left[\begin{array}{c}{S}_{\mathrm{yk}}\hfill \\ \multicolumn{1}{c}{v_{\mathrm{yk}}}\\ \multicolumn{1}{c}{{\phi }_k}\\ \multicolumn{1}{c}{{\stackrel{·}{\phi }}_k}\\ \multicolumn{1}{c}{y_{\mathrm{ref},k}}\end{array}\right]{\mathit{A}}_{\mathit{d}}=\left[\begin{array}{ccccc}1\hfill & T_s\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & T_s\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & I\hfill \end{array}\right]{\mathit{B}}_{\mathit{d}}=\left[\begin{array}{cc}0.5{T_s}^2\hfill & 0\hfill \\ \multicolumn{1}{c}{T_s}& 0\hfill \\ \multicolumn{1}{c}{0}& 0.5{T_s}^2\hfill \\ \multicolumn{1}{c}{0}& T_s\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]$$

(22)

$${\mathit{C}}_{\mathit{d}}=\left[\begin{array}{ccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & g\hfill & 0\hfill & -1\hfill \end{array}\right]{\mathit{D}}_{\mathit{d}}=\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill \end{array}\right]{\mathit{u}}_{\mathit{k}}=\left[\begin{array}{c}{a}_k\hfill \\ \multicolumn{1}{c}{{\stackrel{··}{\phi }}_k}\end{array}\right],$$

(23)

the braking law is defined as

$$\begin{array}{c}|S_{\mathrm{yk}}+c_v{T}_1{\stackrel{·}{S}}_{\mathrm{yk}}+0.5c_u{\stackrel{··}{S}}_{\mathrm{yk}}{T}_1^2|<0.5S_{y,\max }\hfill \\ |{\phi }_k+c_v{T}_2{\stackrel{·}{\phi }}_k+0.5c_u{\stackrel{··}{\phi }}_k{T}_2^2|<0.5{\phi }_{\max },\hfill \end{array}$$

(24)

which subject to equation (25)

$$\begin{array}{c}\left|x_{1k}\right|\le S_{y,\max }\left|x_{2k}\right|\le v_{y,\max }\left|u_{1k}\right|\le a_{y,\max }\hfill \\ \left|x_{3k}\right|\le {\phi }_{\max }\left|x_{4k}\right|\le {\omega }_{\max }\left|u_{2k}\right|\le {\stackrel{·}{\omega }}_{\max },\hfill \end{array}$$

(25)

where ${\phi }_k$ , ${\stackrel{·}{\phi }}_k$ , $v_{\mathrm{yk}}$ , and $S_{\mathrm{yk}}$ , represent for the discrete state variables of tilted angle, and tilted rate, velocity, and position of the driving simulator, respectively. The reference vector ${\mathit{y}}_{\mathit{ref}{,}_{\mathit{k}}}$ (Figure 8) is selected for pulling the simulator to its home position while tracking the target specific force $f_{\mathrm{vy}}$ . Similarly, the discrete acceleration $a_{\mathrm{sy}}$ and angular acceleration ${\stackrel{··}{\phi }}_{\mathrm{sx}}$ are $a_k$ and ${\stackrel{··}{\phi }}_k$ , respectively. The limits values $\{.{\}}_{\max }$ , are defined as the absolute maximum values of the simulator’s motion quantities such as position, velocity, lateral acceleration, and the motion perception threshold: angular velocity, and angular acceleration. The sample time of the discrete dynamic model is $T_s$ , and the parameters of the braking control law including $c_v$ , $c_u$ , $T_1$ , and $T_2$ determine the smoothing approach of the simulator to its workspace limits. Additionally, the important condition for selecting possible parameters in order to avoid overshoots in equation (24) is recommend as $(c_v^2/2c_u)>1$ . Starting from the original parameters published by Fang and Kemeny, ¹² the values of parameters are selected to gain feasible simulations for tracking specific force on the Robocoaster, obtaining the “exMPC” algorithm (Table 10).

Figure 8.

Control diagram of the explicit MPC.

Table 10.

Explicit MPC parameters as tuned for the KUKA robocoaster motion platform.

Model		Limit values		Parameters of law brake
$T_s$	0.125 (s)	$S_{y\max }$	1.5 (m)	$c_u$	1.5 (s)
$q_{\mathrm{yi}}$	$4|0.01|299$	$v_{y\max }$	10 (m/s)	$c_v$	1 (s)
$r_i$	$0.1|2$	$a_{y\max }$	14 (m/s2)	$T_1$	2 (s)
$q_{\mathrm{xi}}$	$1|1|1|1$	${\phi }_{\max }$	0.5 (rad)	$T_2$	1 (s)
$N$	5	${\omega }_{\max }$	0.25 (rad/s)	–	–
$\mathrm{Nc}$	5	${\stackrel{·}{\omega }}_{\max }$	10 (rad/s2)	–	–

Numerical comparison of MCA families

In order to compare the perceived quality of outcomes of the different MCA families, a simple virtual ride at a constant speed of $6$ m/s along a planar left-right S-shaped trajectory with smoothed transitions was chosen (thick lines in Figure 9), featuring purely lateral accelerations. Hereby, two cases were considered. In case 1, the original target lateral acceleration was fed into the $11$ reviewed MCAs. In case 2, for a preparatory subjective test, which verifies the correlation between the numerical criteria, and passengers’ rates according to their motion perception, the target accelerations were scaled according to the rules a) the range of simulated motion is within the workspace of the simulators, b) the $k_y=[0.4,1]$ (equation (27)). The range of $k_y$ was referred to the scale factor experiment of Berthoz et al. ²³ The algorithm MPC was removed from this case as it featured a too large phase shift with respect to the original signal in the specific force.

Figure 9.

Results of MCA families of lateral specific forces and angular velocity for the test track – case 1: (a) simulated lateral specific forces and (b) simulated angular velocites.

Figure 9(a) for case 1 shows that most signals follow well the target specific force, except the original MPC, which displays a strong phase shift. On the other hand, Figure 9(b) shows that all MCAs violate the threshold level for the angular velocity, beyond which the tilt rotation is sensed as such and not as an increase of lateral g-force. Hereby, the most severe angular velocity threshold violation is reached by the exMPC algorithm. By scaling down the target accelerations (Figure 10), the angular velocity violations could be reduced for all MCAs except the exMPC algorithm, but the target specific force was not reached by all MCAs except MPC*.

Figure 10.

Results of MCA families of lateral specific forces and angular velocity for the test track – case 2: (a) simulated lateral specific forces and (b) simulated angular velocites.

For the motion perception evaluation, only a few objective measures exist in the literature. A comprehensive objective measures were proposed by the “performance indicator” of Pouliot et al. ¹⁹ and the “good criterion” of Fischer. ²⁰

The “performance indicator,” ${\lambda }_1$ and ${\lambda }_2$ (equation (26)), analyzes the quality of the simulated signals. The first indicator ${\lambda }_1$ is the normalized average of the square error between the target and simulated signals, while the second indicator ${\lambda }_2$ is the normalized average square error of the rates of change. Each indicator combines two sub-indicator indexed “ $f$ ” for the specific force, and indexed “ $\omega$ ” for the angular velocity. The two sub-indicators are normalized by $a_{\max }$ and ${\omega }_{\max }$ (equation (26)) which represent the maximal acceleration and the maximal angular velocity of the simulator, respectively.

$${\lambda }_1=100·(\frac{{\lambda }_{1,f}}{a_{\max }}+\frac{{\lambda }_{1,\omega }}{{\omega }_{\max }});{\lambda }_2=100·(\frac{{\lambda }_{2,f}}{a_{\max }}+\frac{{\lambda }_{2,\omega }}{{\omega }_{\max }})$$

(26)

Further, the “good criterion,” ²⁰ separates the translational and rotational sub-indicators into a “scale” component (index “sc”) for estimating strong stimulation of simulated signals, and a “shape” component (index “sh”) for predicting the false cues in the simulated signals. We describe here the scale and shape components for the specific force that applies analogously to the rotational part. Considering on the time $t_j$ in the N + 1 time sample, let $f_{\mathrm{si}}^j,f_{\mathrm{vi}}^j$ , be the simulated and target (“vehicle”) specific forces, respectively, for a the coordinate $i\in \{x,y,z\}$ . Regarding the simulated and target specific force, the global scaling factor is defined as equation (27)

$$k_i=(\sum _{j=0}^N|f_{\mathrm{si}}^j|)/(\sum _{j=0}^N|f_{\mathrm{vi}}^j|),i\in \{x,y,z\}.$$

(27)

At a sample point, the different components are equation (28).

$$\mathrm{\Delta }{f}_{i,\mathrm{sc}}^j=f_{\mathrm{vi}}^j(1-k_i);\mathrm{\Delta }{f}_{i,\mathrm{sh}}^j=k_i{f}_{\mathrm{vi}}^j-f_{\mathrm{si}}^j,i\in \{x,y,z\}$$

(28)

The partial criteria for specific force components

$${\lambda }_{1f,\mathrm{sc}}=\frac{1}{N+1}·\sum _{j=0}^N\sqrt{\sum _{i\in \{x,y,z\}}{(\mathrm{\Delta }{f}_{i,\mathrm{sc}}^j)}^2}$$

(29)

$${\lambda }_{1f,\mathrm{sh}}=\frac{1}{N+1}·\sum _{j=0}^N\sqrt{\sum _{i\in \{x,y,z\}}{(\mathrm{\Delta }{f}_{i,\mathrm{sh}}^j)}^2}$$

(30)

From this, the “good criteria” for simulator cues and their time derivative are obtained as:

$${\lambda }_1^{*}=(100/f_{\mathrm{nr}})·{\lambda }_{1f}+(100/{\omega }_{\mathrm{nr}})·{\lambda }_{1\omega }$$

(31)

$${\lambda }_2^{*}=(100/{\stackrel{·}{f}}_{\mathrm{nr}})·{\lambda }_{2f}+(100/{\stackrel{·}{\omega }}_{\mathrm{nr}})·{\lambda }_{2\omega }$$

(32)

With

$${\lambda }_{\mathrm{pq}}=({\lambda }_{\mathrm{pq},\mathrm{sc}}+{\lambda }_{\mathrm{pq},\mathrm{sh}}),p\in \{1,2\};q\in \{f,\omega \}$$

(33)

$$f_{\mathrm{nr}}=\max (|f_{\mathrm{vi}}|)(i\in \{x,y,z\})$$

$${\omega }_{\mathrm{nr}}=\max (|{\omega }_{\mathrm{vx}}+{\omega }_{\mathrm{sx}}|,|{\omega }_{\mathrm{vy}}+{\omega }_{\mathrm{sy}}|),|{\omega }_{\mathrm{sz}}|)$$

$${\stackrel{·}{f}}_{\mathrm{nr}}=\max (|{\stackrel{·}{f}}_{\mathrm{vi}}|)(i\in \{x,y,z\})$$

$${\stackrel{·}{\omega }}_{\mathrm{nr}}=\max (|{\stackrel{·}{\omega }}_{\mathrm{sx}}+{\stackrel{·}{\omega }}_{\mathrm{vx}}|,|{\stackrel{·}{\omega }}_{\mathrm{sy}}+{\stackrel{·}{\omega }}_{\mathrm{vy}}),|{\stackrel{·}{\omega }}_{\mathrm{sz}}|)$$

The comparison the two objective measures are shown in Figure 11(a), and the results of the partial “good” components in Figure 11(b), both for case 2, as case 1 produced too high false angular velocity cues for almost all rides.

Figure 11.

Comparison of individual good criterion metrics – case 2: (a) performance indictors and (b) good criterion.

One can see that both objective measures yield a similar distribution for the different MCAs, and that the path-tuned optimal control ZyRo* algorithm displays the second-best perception quality. The algorithm has better both index of the good criterion than explicit model-control exMPC algorithm, which take computation times of approximately from 30 min to 3 h on a PC desktop with CPU type—quadcore 2.66 GHz, as compared to 1.5 min for the ZyRo* algorithm. The poor performance of the exMPC algorithm may be explained by the required handling of the ride simulator limits, which forces the simulator to decelerate even when the tilting angle is not yet large enough to compensate for the lack of specific force, leading to unrealistic departures from the target acceleration. Overall the here newly proposed adapted MPC* and ZyRo* algorithms are shown to yield the best results for both objective measures. However, when the computational complexity of the MPC* and the ZyRo* are compared, the MPC* consumes roughly thirty times the ZyRo*. The length of the prediction horizon and the complexity of the control model used in the MPC* appear to be the reasons behind the large time comsuming. Moreover, the parameters of the ZyRo algorithm are transparent to the simulated quantities, making parameter adjustment easier than that of classical, adaptive, and optimum washout filters, as well as MPC*. However, because the washout effect is not used in the simulated angles, the steady-state must be adjusted to bring the final simulated angle to zero.

It’s worth noting that both MPC* and ZyRo* effectively use the motion platform’s workspace to discover the best combination of translational and rotational motion to reproduce the two biggest scaled specific forces. However, because the MPC* closely matches the target specific force in Figure 10(a), it provides angular velocity false cues.

Preliminarily subjective validation

The objective criteria normally provide the first estimation of the simulated motion fidelity regarding the predicted false cues. Considering the human perception, the second phase test of simulated motion with the real passengers was implemented in the next-step on the motion platform. The preliminary subjective tests are important to validate the correlation between objective and subjective measures. Similar tests were performed in Berthoz et al. ²³ to investigate the effect of the washout-filters. The test was implemented with a total of 17 participants including 12 males and 5 females who have the age in range 20–30 years. The participants were asked to rate the 10 rides of case 2 on five satisfying levels in the range {very good, good, satisfactory, sufficient, not sufficient} corresponding to the grade scale $\{1,2,3,4,5\}$ with $\pm 0.3$ as shown in the Table 11.

Table 11.

Subjective grade scale.

Grade	1	1.3	1.7	2	2.3	2.7	3	3.3	3.7	4	4.3	4.7	5
Meaning	Very good		Good			Satisfactory			Sufficient		Not sufficient

This scaling based on German school grading marks was used due to well-perceived by all test subjects and no more guide was needed. A participant successively tried 10 rides in arbitrary order and asked to give his/her evaluation after each ride. The interclass correlation for rates of 17 judges is ICC(3,k) = 0.4937 with p = 0.0462—probability of false-rejection—(k = 17 according to Berthoz et al. ²³ ). The ICC values express the fair reliability of the rates.

The individual subjective ratings (SS) and the scale factor (SF) of simulated specific forces are shown in Table 12. For the graphical representation, the “bubble plot” (Figure 12) uses the circle size to describe the number of rates for individual subjective scores which are fairly variable in the full scale 1–5, in general. However, there are also some abnormal rates for such worst ride (exMPC algorithm) that one participant gives the best score as “very good.” The abnormal scores are controversial and could be happened because the participants have the particular threshold of rotational motion cue.

Table 12.

Average subjective scores.

	CL	AD	Tel-YoMe	Tel	OpS	OpRN	MPC*	exMPC	ZyRo	ZyRo*
SS	2.24	2.35	2.24	2.12	2.06	2.06	2.24	2.82	2.47	2.18
SF	0.5	0.5	0.42	0.42	0.42	0.42’	1	0.58	0.4	0.58

Figure 12.

Individual subjective and objective score.

In order to compare the objective perception ratings with the subjective ratings, the partial components of the objective ratings were fitted by linear regression with the average subjective ratings. Besides, we introduced as an additional component for the objective rating a penalty-function $J_{\omega }$ describing the false cue due to the violation of the angular velocity threshold as

$$J_{\omega }=\{\begin{array}{cc}0\hfill & \mathrm{if}\max (\left|{\omega }_s\right|)<{\delta }_S\hfill \\ \multicolumn{1}{c}{e^{{(\max (\left|{\omega }_s\right|)-{\delta }_S)}^2}-1}& \mathrm{if}\max (\left|{\omega }_s\right|)\ge {\delta }_S\hfill \end{array}$$

(34)

In a first round, we fitted the parameters ${\lambda }_1$ and ${\lambda }_2$ of both objective measures together with the penalty function $J_{\omega }$ to the average subjective measures for the linear regression

$$F=w_0+w_{{\lambda }_1^{*}}{\lambda }_1^{*}+w_{{\lambda }_2^{*}}{\lambda }_2^{*}+w_{J_{\omega }}{J}_{\omega }$$

(35)

with the regress function in MATLAB. In order to have an idea of the impact of each term, we switched on and off each term for all possible combinations. The results are shown in the Table 13, which displays the root mean square error (RMSE) and Pearson correlation (CORR) for each combination of terms. One can see that both ${\lambda }_1$ and ${\lambda }_2$ have a bad correlation with the average subjective measure and that even the best fit of ${\lambda }_1$ nor ${\lambda }_2$ together yield only a fair correlation. However, by adding the penalty term $J_{\omega }$ , this term alone gives a better correlation with the average subjective measure than the best of ${\lambda }_1$ and ${\lambda }_2$ , which is improved only slightly if the aforementioned two are added to the linear regression. This shows that false angular velocity cues have a significant impact on motion perception quality.

Table 13.

Coefficient of numerical expected functions $\{F_1,\cdots ,F_7\}$ .

	$w_0$	$w_{{\lambda }_1^{*}}$	$w_{{\lambda }_2^{*}}$	$w_{J_{\omega }}$	RMSE	CORR
$F_1$	2.2716	0.0075	0	0	0.2167	0.0056
$F_2$	2.4314	0	−0.2199	0	0.2091	0.2624
$F_3$	1.8270	1.2780	−0.9127	0	0.1825	0.5396
$F_4$	2.2076	0	0	0.6145	0.1184	0.8378
$F_5$	2.1985	0.0107	0	0.6145	0.1183	0.8378
$F_6$	2.2100	0	−0.0032	0.6136	0.1183	0.8378
$F_7$	2.1722	0.0910	−0.0577	0.5990	0.1182	0.8382

In a second round, we split the linear regression model even finer by taking into account all partial components of the “good criterion” as

$$F=\sum _{q\in \{1,2\}}\sum _{p\in \{f,\omega \}}\sum _{n\in \{\mathrm{sc},\mathrm{sh}\}}{w}_{\mathrm{qp},n}·{\lambda }_{\mathrm{qp},n}+w_{J_{\omega }}{J}_{\omega }$$

(36)

The coefficients $w_{\mathrm{qp},n}$ are listed in the Table 14. One can see that applying the suitable fitting coefficients to the “Good criterion” components leads very good correlations with the average subjective score achieved. This is even further optimized when including the penalty term $J_{\omega }$ , yielding Pearson correlation $95.42\mathrm{\%}$ and a root mean square error of 0.0660. By plotting the thus optimized numerical scores in the individual/average subjective rating plot of Figure 12 (red crosses). It could be noted that the numerical scores track much better the average subjective score than the individual ones. This is seen also in Figure 13, which shows the ranking of the 10 MCA both in terms of average subjective and objective measures.

Table 14.

Coefficient of numerical expected functions $\{F_8,F_9\}$ .

	$w_{1f,\mathrm{sc}}$	$w_{2f,\mathrm{sc}}$	$w_{1f,\mathrm{sh}}$	$w_{2f,\mathrm{sh}}$	$w_{1\omega ,\mathrm{sh}}$	$w_{2\omega ,\mathrm{sh}}$	$w_{J_{\omega }}$	RMSE	CORR
$F_8$	14.6490	−13.7195	−0.9510	2.6956	0.1424	−0.2323	0	0.0788	0.9352
$F_9$	9.1670	−7.9579	−0.2280	1.8073	1.9510	−2.1366	1.4126	0.0660	0.9542

Figure 13.

MCA ranking by average subjective and numerical score.

The preliminary results thus show that (1) the individual objective parameters ${\lambda }_1$ and ${\lambda }_2$ seem not to correlate well with average subjective ratings; (2) the violation of the false angular velocity cue threshold influences significantly the subjective ride perception; (3) after down-scaling target accelerations there are little perception differences for most MCAs; and (4) by proper linear regression of the partial components of the “good criterion” together with the penalty terms for false cue angular velocity a very good correlation between average subjective and objective measures can be achieved. However, further research is necessary to better understand the factors impacting on subjective ratings of a simulator task because the individual scores are not consistent and have fair reliability.

Conclusions

The present paper first reviews the families of offline MCAs from the literature including the classical washout, adaptive washout, optimal washout, optimal control, and model-predictive control MCAs. Then, the MCAs were implemented in MATLAB for virtual rides over a planar S-shaped curve with lateral acceleration only. Next, the comparison of their behavior and the performance evaluation are carried out using both “performance indicator” and “good criterion” in which the deviation of the simulated and target motion cues (specific forces and angular velocities) are estimated. The comparison of different simulated responses in the case $1$ shows that the optimal tracking control approach (ZyRo algorithm) with the tuning parameters yields the best results in tracking the target signals. The ZyRo and MPC algorithm exploit better workspace of the simulator. In the comparison of case $2$ and the preliminary subjective evaluation, because of the clear disagreement of current numerical criteria and subjective evaluation, the proposed numerical criteria are built as the linearized function of performance indexes or sub-indexes. The penalty function $J_{\omega }$ , inserted in some functions to increase the correlation of numerical score w.r.t subjective score, shows the significant effect of false cues of angular velocity to the motion fidelity. The function $F_9$ has the best correlation and could be used as the numerical measure for the numerical prediction of sensations. In this paper, only the numerical measures of the simulated cues, that are generated by motion cueing algorithms, and the preliminary results with 17 test-subjects have been addressed. The proposed numerical criteria can also be utilized to evaluate the MCA response’s performance quality throughout the parameter tuning process, reducing the amount of time spent on trial and error.