A Human Motion Compensation Framework for
a Supernumerary Robotic Arm

Here we briefly introduce the general layout of our SRA. The detailed description is provided in our previous work [10]. The 6-DoF robotic arm is mounted on the user’s shoulder, which aims to guarantee the maximum collision-free workspace between the SRA and the human arm. The unique advantage of our SRA is the body adaptability. We design an adjustable mechanical interface allowing different users to manually adjust the installation position of the manipulator to fit different users’ body dimensions and switch left/right shoulders. The soft shoulder support keeps the SRA stable and comfortable on the user’s body.

To achieve the kinematic chain integration of the SRA and the human, the IMU-based motion capture system (Xsens MVN) is chosen as the sensory interface due to its portability. This motion capture system combines the biomechanical model of the human user and an IMU sensor network to estimate the 3D kinematic data of the user. First, the user body dimensions should be measured, including body height, foot size, arm span, ankle height, hip height, hip width, knee height, shoulder width, and shoe sole height. Then, the user performs a calibration process wearing a flexible suit and body straps to attach 17 IMU wireless tracker modules, and the IMU signal data is transmitted to the master (the transmission station and the USB dongle). Through our Xsens-to-ROS bridge [11], the user’s motion is detected in real time, which works as an important feedback for controlling the SRA.

The velocity-level kinematics of the floating-base system[12] can be established as follows:

$\displaystyle\boldsymbol{\dot{x}}_{E}=\boldsymbol{J}_{B}\boldsymbol{\dot{x}}_{B}+\boldsymbol{J}_{M}\boldsymbol{\dot{\theta}_{M}}$

(1)

where $\boldsymbol{\dot{x}}_{E}=\left[\boldsymbol{v}_{E}^{T},\boldsymbol{\omega}_{E}^{T}\right]^{T}\in\mathbb{R}^{6}$ is the generalized velocity of the EE, $\boldsymbol{\dot{x}}_{B}=\left[\boldsymbol{v}_{B}^{T},\boldsymbol{\omega}_{B}^{T}\right]^{T}\in\mathbb{R}^{6}$ is the generalized velocity of the human base, $\boldsymbol{\theta}_{M}=\left[\theta_{1},\theta_{2},\cdots,\theta_{6}\right]^{T}\in\mathbb{R}^{6}$ is the joint position of the manipulator, and $\boldsymbol{J}_{B}\in\mathbb{R}^{6\times 6}$ and $\boldsymbol{J}_{M}\in\mathbb{R}^{6\times 6}$ are Jacobian matrices of the base and the manipulator defined in $\Sigma_{I}$. $\boldsymbol{J}_{B}$ is defined as

$$\boldsymbol{J}_{B}=\begin{bmatrix}\boldsymbol{E}_{3\times 3}&-\boldsymbol{\widetilde{p}}_{BE}\\ \boldsymbol{0}_{3\times 3}&\boldsymbol{E}_{3\times 3}\end{bmatrix}$$

(2)

where $\boldsymbol{E}$ is the identity matrix. The symbol $\left(\tilde{\cdot}\right)$ denotes a skew-symmetric matrix, namely, for an arbitrary vector $\boldsymbol{k}=\left[k_{x},k_{y},k_{z}\right]^{T}\in\mathbb{R}^{3}$, $\tilde{\boldsymbol{k}}$ is defined as

$\displaystyle\tilde{\boldsymbol{k}}=\begin{bmatrix}0&-k_{z}&k_{y}\\ k_{z}&0&-k_{x}\\ -k_{y}&k_{x}&0\end{bmatrix}$

(3)

and $\boldsymbol{p}_{BE}\in\mathbb{R}^{3}$ is the position vector from $\Sigma_{B}$ to $\Sigma_{E}$ in $\Sigma_{I}$:

$$\boldsymbol{p}_{BE}=\boldsymbol{R}_{B}\cdot\boldsymbol{p}_{BE}^{B}$$

(4)

where $\boldsymbol{R}_{B}\in\mathbb{R}^{3}$ is the orientation of $\Sigma_{B}$ w.r.t. $\Sigma_{I}$, and $\boldsymbol{p}_{BE}^{B}\in\mathbb{R}^{3}$ is the position vector from $\Sigma_{B}$ to $\Sigma_{E}$ in $\Sigma_{B}$. Meanwhile, we can deduce

$$\boldsymbol{J}_{M}=\begin{bmatrix}\boldsymbol{R}_{B}&\boldsymbol{0}_{3\times 3}\\ \boldsymbol{0}_{3\times 3}&\boldsymbol{R}_{B}\end{bmatrix}\cdot\boldsymbol{J}_{M}^{B}$$

(5)

where $\boldsymbol{J}_{M}^{B}\in\mathbb{R}^{6\times 6}$ is the Jacobian matrix of the manipulator defined in $\Sigma_{B}$. According to the assumed precondition (P1), we know $\boldsymbol{\omega}_{B}\approx\boldsymbol{0}$ and $\boldsymbol{J}_{B}\boldsymbol{\dot{x}}_{B}=[\boldsymbol{v}_{B}^{T},\mathbf{0}^{T}]^{T}$ when the human user is performing the translational motion at the local working position. Therefore, we can obtain a simplified kinematics model:

$$\begin{bmatrix}\boldsymbol{v}_{E}\\ \boldsymbol{\omega}_{E}\end{bmatrix}=\begin{bmatrix}\boldsymbol{v}_{B}\\ \boldsymbol{0}\end{bmatrix}+\begin{bmatrix}\boldsymbol{J}_{Mv}\\ \boldsymbol{J}_{M\omega}\end{bmatrix}\boldsymbol{\dot{\theta}_{M}}$$

(6)

where $\boldsymbol{J}_{Mv}$ and $\boldsymbol{J}_{M\omega}\in\mathbb{R}^{3\times 6}$ are the translational and rotational Jacobian matrices of $\boldsymbol{J}_{M}$, respectively.

Our basic objective is to keep the EE stable when the user is moving. Due to the translational motion of the human ($\boldsymbol{v}_{B}\neq\boldsymbol{0}$), we compensate for $\boldsymbol{v}_{B}$ by planning $\boldsymbol{\dot{\theta}_{M}}$ to keep $\boldsymbol{x}_{E}$ still (or to follow a desired trajectory). As our manipulator does not have redundant DoFs, it’s hard to keep the 6D pose of the EE unchanged at any point in the whole workspace. Since this problem cannot be addressed, we make a compromise by prioritizing subtasks where we assign the EE position compensation ($\boldsymbol{v}_{E}=0$) as the primary task and the EE orientation compensation ($\boldsymbol{\omega}_{E}=0$) as the secondary task. According to (6), the primary task of the compensation problem is formulated as

$\displaystyle\underset{\boldsymbol{\dot{\theta}}_{M}}{\text{min}}~{}~{}~{}$

$\displaystyle\|\boldsymbol{J}_{Mv}\boldsymbol{\dot{\theta}}_{M}+\boldsymbol{v}_{B}\|^{2}$

(7)

Moreover, we assume $\boldsymbol{\dot{\theta}}_{M}^{\star}\in\mathbb{R}^{6}$ is the solution of (7) and the secondary task is to meet

	$\displaystyle\underset{\boldsymbol{\dot{\theta}}_{M}}{\text{min}}~{}~{}~{}$	$\displaystyle\|\boldsymbol{J}_{M\omega}\boldsymbol{\dot{\theta}}_{M}\|^{2}$		(8a)
	s.t	$\displaystyle\boldsymbol{J}_{Mv}\boldsymbol{\dot{\theta}}_{M}=\boldsymbol{J}_{Mv}\boldsymbol{\dot{\theta}}_{M}^{\star}$		(8b)

which means the solution of (8) is within the set of optimizing the primary task [13]. Meanwhile, the planned trajectories should meet the physical constraints ($\boldsymbol{\theta}_{min}\leq\boldsymbol{\theta}_{M}\leq\boldsymbol{\theta}_{max}$ and $\boldsymbol{\dot{\theta}}_{min}\leq\boldsymbol{\dot{\theta}}_{M}\leq\boldsymbol{\dot{\theta}}_{max}$), where $\boldsymbol{\theta}_{min}$ and $\boldsymbol{\theta}_{max}\in\mathbb{R}^{6}$ are the lower and upper boundaries of the joint position vector; $\boldsymbol{\dot{\theta}}_{min}$ and $\boldsymbol{\dot{\theta}}_{max}\in\mathbb{R}^{6}$ are the lower and upper boundaries of the joint velocity vector.

In this subsection, two approaches are used to solve the above problem: 1) the nullspace-based method (NBM) and 2) the reconstructed Jacobian method (RJM).

The human motion compensation framework is given in Fig. 3, which includes 5 parts: (a) the physical integration, (b) the kinematic chain integration, (c) the coordinate frame transformation, (d) the motion planning, and (e) the motion control. First, the human user wears the SRA with the IMU-based sensory interface to establish the human-robot integration system. After the calibration, the inertial frame $\Sigma_{I}$ lies in the right foot of the user at the initial calibration position (Fig. 2) and all the kinematic information of the human user is defined w.r.t. $\Sigma_{I}$. As the user usually works at a local working position, the compensation motion control should be finished at the local position. Therefore, we set a local home frame and the human motion should be defined w.r.t. the local home frame.

As shown in Fig. 3b, we aim to achieve the motion compensation at a local working position, the human torso frame $\Sigma_{H}$ at $t_{0}$ (noted as $\Sigma_{H,t_{0}}$) is selected as the local home frame. Therefore, the human translation motion should be defined w.r.t. $\Sigma_{H,t_{0}}$. $\Sigma_{H,t_{0}}$ is set as follows: First, the human user reaches the local working position and stands firmly. Then, we record the local torso frame $\Sigma_{H,t_{0}}$ with its current position $\boldsymbol{p}_{H,t_{0}}$ and its current orientation $\boldsymbol{R}_{H,t_{0}}^{I}$ w.r.t. $\Sigma_{I}$, respectively. According to the assumed precondition (P1), we can deduce the position variation of the local home frame $\Sigma_{H}$ from $t_{0}$ to $t_{k}$ w.r.t. $\Sigma_{H,t_{0}}$ as

$$\Delta\boldsymbol{p}{{}_{H,t_{k}}^{H,t_{0}}}=\boldsymbol{R}_{I}^{H,t_{0}}(\boldsymbol{p}_{H,t_{0}}-\boldsymbol{p}_{H,t_{k}})$$

(16)

where $\boldsymbol{p}_{H,t_{k}}$ is the position of $\Sigma_{H}$ at $t_{k}$ w.r.t. $\Sigma_{I}$. Similarly, we can obtain the local position variation of the base frame $\Sigma_{B}$ from $t_{0}$ to $t_{k}$ w.r.t. $\Sigma_{H,t_{0}}$ as

$$\Delta\boldsymbol{p}{{}_{B,t_{k}}^{H,t_{0}}}=\boldsymbol{R}_{I}^{H,t_{0}}(\boldsymbol{p}_{H,t_{0}}+\boldsymbol{p}_{HB,t_{0}}-\boldsymbol{p}_{H,t_{k}}-\boldsymbol{p}_{HB,t_{k}})$$

(17)

where $\boldsymbol{p}_{HB,t_{0}}$ and $\boldsymbol{p}_{HB,t_{k}}\in\mathbb{R}^{3}$ are the relative position vectors from $\Sigma_{H}$ to $\Sigma_{B}$ at $t_{0}$ and $t_{k}$ w.r.t. $\Sigma_{I}$, respectively. According to the assumed precondition (P2), we know $\boldsymbol{p}_{HB,t_{0}}=\boldsymbol{p}_{HB,t_{k}}$, so we can deduce

$$\Delta\boldsymbol{p}{{}_{B,t_{k}}^{H,t_{0}}}=\Delta\boldsymbol{p}{{}_{H,t_{k}}^{H,t_{0}}}$$

(18)

$$\boldsymbol{v}{{}_{B,t_{k}}^{H,t_{0}}}=\boldsymbol{v}{{}_{H,t_{k}}^{H,t_{0}}}=\boldsymbol{R}_{I}^{H,t_{0}}\boldsymbol{v}_{H,t_{k}}$$

(19)

where $\boldsymbol{v}_{H,t_{k}}\in\mathbb{R}^{3}$ is the velocity vector of $\Sigma_{H}$ at $t_{k}$ w.r.t. $\Sigma_{I}$, $\boldsymbol{v}{{}_{B,t_{k}}^{H,t_{0}}}$ and $\boldsymbol{v}{{}_{H,t_{k}}^{H,t_{0}}}\in\mathbb{R}^{3}$ are velocity vectors of $\Sigma_{B}$ and $\Sigma_{H}$ at $t_{k}$ w.r.t. $\Sigma_{H,t_{0}}$, respectively. Moreover, we can deduce the local position variation of $\Sigma_{E}$, namely $\Delta{\boldsymbol{p}_{E}}$ at $t_{k}$ w.r.t. $\Sigma_{H,t_{0}}$, which will be used to control our SRA:

$$\Delta\boldsymbol{p}{{}_{E,t_{k}}^{H,t_{0}}}=\Delta\boldsymbol{p}{{}_{B,t_{k}}^{H,t_{0}}}+\Delta\boldsymbol{p}{{}_{BE,t_{k}}^{H,t_{0}}}$$

(20)

where $\Delta\boldsymbol{p}{{}_{BE,t_{k}}^{H,t_{0}}\in\mathbb{R}^{3}}$ is the relative position from $\Sigma_{B}$ to $\Sigma_{E}$ at $t_{k}$ w.r.t. $\Sigma_{H,t_{0}}$ and expressed as

$$\Delta\boldsymbol{p}_{{BE,t_{k}}}^{H,t_{0}}=\boldsymbol{R}_{H,t_{k}}^{H,t_{0}}\cdot\boldsymbol{R}_{B,t_{k}}^{H,t_{k}}\cdot\boldsymbol{p}_{E,t_{k}}^{B,t_{k}}-\boldsymbol{R}_{B,t_{0}}^{H,t_{0}}\cdot\boldsymbol{p}_{E,t_{0}}^{B,t_{0}}$$

(21)

where $\boldsymbol{R}_{H,t_{k}}^{H,t_{0}},\boldsymbol{R}_{B,t_{k}}^{H,t_{k}}$, and $\boldsymbol{R}_{B,t_{0}}^{H,t_{0}}\in\mathbb{R}^{3\times 3}$ are the rotation matrices of $\Sigma_{H,t_{k}}$ w.r.t. $\Sigma_{H,t_{0}}$, $\Sigma_{B,t_{k}}$ w.r.t. $\Sigma_{H,t_{k}}$, and $\Sigma_{B,t_{0}}$ w.r.t. $\Sigma_{H,t_{0}}$, respectively. According to the assumed preconditions (P1) and (P2), each of them is $\boldsymbol{E}_{3\times 3}$ and (21) can be reduced into

$$\Delta\boldsymbol{p}_{{BE,t_{k}}}^{H,t_{0}}=\boldsymbol{p}_{E,t_{k}}^{B,t_{k}}-\boldsymbol{p}_{E,t_{0}}^{B,t_{0}}$$

(22)

For the EE orientation stabilization, controlling it w.r.t. $\Sigma_{B}$ is the same as controlling it w.r.t. $\Sigma_{H}$, so we use the following quaternion error to control our SRA:

$$\Delta\boldsymbol{Q}_{E,t_{k}}^{B}=\boldsymbol{Q}_{E,t_{0}}^{B,t_{0}}\otimes\boldsymbol{Q}_{E,t_{k}}^{B,t_{k}-1}=\left\{\Delta\eta_{E,t_{k}}^{B},\Delta\boldsymbol{\epsilon}_{E,t_{k}}^{B}\right\}$$

(23)

When the manipulator is close to a singularity configuration, the inverse kinematics becomes unstable and the joint increment becomes large. The singular value filtering (SVF) method [16] is adopted in our control framework, which can guarantee the lower-bounded singular values. According to the singular value decomposition, we can deduce

$$\boldsymbol{J}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{T}=\sum_{i=1}^{r}\sigma_{i}\boldsymbol{u}_{i}\boldsymbol{v}_{i}^{T}$$

(24)

where $\boldsymbol{v}_{i}$ and $\boldsymbol{u}_{i}$ are the input and output singular vectors which are the $i$th columns of $\boldsymbol{U}\in\mathbb{R}^{m\times m}$ and $\boldsymbol{V}\in\mathbb{R}^{n\times n}$, respectively, $\sigma_{i}$ is the singular value meeting $\sigma_{1}\geq\sigma_{2}\geq...\geq\sigma_{r}>0$, with $r\leq$ min$\{m,n\}$ being the rank of $\boldsymbol{J}$. The pseudo-inverse of $\boldsymbol{J}$ is $\boldsymbol{J}^{\dagger}=\boldsymbol{U}\boldsymbol{\Sigma}^{\dagger}\boldsymbol{V}^{T}=\sum_{i=1}^{r}\frac{1}{\sigma_{i}}\boldsymbol{u}_{i}\boldsymbol{v}_{i}^{T}$. The general damped least-squares pseudo-inverse is expressed as

$$\boldsymbol{J}^{\dagger}_{D}=\sum_{i=1}^{r}\frac{\sigma_{i}}{\sigma_{i}^{2}+\lambda^{2}}\boldsymbol{u}_{i}\boldsymbol{v}_{i}^{T}$$

(25)

where $\lambda$ is the damping factor, but a constant small value for $\lambda$ cannot always guarantee good performance in the workspace. Using the SVF method, the filtered Jacobian matrix is defined as

	$$\displaystyle\boldsymbol{J}_{F}=\sum_{i=1}^{r}f(\sigma_{i})\boldsymbol{u}_{i}\boldsymbol{v}_{i}^{T}$$		(26)
	$$\displaystyle f(\sigma)=\frac{\sigma^{3}+\upsilon\sigma^{2}+2\sigma+2\sigma_{0}}{\sigma^{2}+\upsilon\sigma+2}$$		(27)

where $f(\sigma)$ is the filtering function, $\sigma_{0}$ is the minimum value that we want to impose on the singular values of $\boldsymbol{J}$, $\upsilon$ is a shape factor, and the corresponding pseudo-inverse of $\boldsymbol{J}_{F}$ is

$$\boldsymbol{J}^{\dagger}_{F}=\sum_{i=1}^{r}\frac{1}{f(\sigma_{i})}\boldsymbol{u}_{i}\boldsymbol{v}_{i}^{T}$$

(28)

To avoid joint velocity and angle limits, we use the joint clamping algorithm [17] to constrain the joint increment.

$$\displaystyle\dot{\theta}_{i}=\begin{cases}-\dot{\theta}_{iLimit},&\dot{\theta}_{i}<-\dot{\theta}_{iLimit}\\ \dot{\theta}_{iLimit},&\dot{\theta}_{i}>+\dot{\theta}_{iLimit}\\ \dot{\theta}_{i},&\text{otherwise.}\\ \end{cases}$$

(29)

where $\dot{\theta}_{iLimit}$ is the $i$th joint velocity boundary. For the joint position at step $k$, we know $\boldsymbol{{\theta}}_{M}(k+1)=\boldsymbol{{\theta}}_{M}(k)+\boldsymbol{H}\Delta\boldsymbol{{\theta}}_{M}$, where $\Delta\boldsymbol{{\theta}}_{M}=\boldsymbol{\dot{\theta}}_{M}\Delta t$ and $\boldsymbol{H}=diag(H_{1},H_{2},...,H_{6})\in\mathbb{R}^{6\times 6}$ is a diagonal activation matrix,

$$\displaystyle{H}_{i}=\begin{cases}1,&{\theta}_{iMin}<{\theta}_{i}<{\theta}_{iMax}\\ 0,&\text{otherwise.}\\ \end{cases}$$

(30)

where ${\theta}_{iMin}$ and ${\theta}_{iMax}$ are the $i$th joint position minimum and maximum boundaries, respectively. Through the joint clamping algorithm, we can guarantee that the planned joint trajectories do not violate the joint limits.

We conduct three 1D motion experiments, namely, the up-down (U-D) motion (Fig. 4), the left-right (L-R) motion (Fig. 5), and the forward-backward (F-B) motion (Fig. 6). Although the user tries to achieve the 1D motion, the human body’s movements are inevitably coupled sometimes. Note that we focus on the z-axis stabilization in the U-D motion, the y-axis stabilization in the L-R motion, and the x-axis stabilization in the F-B motion. As shown in Fig. 7, the human user moves randomly in the 3D motion experiments and we focus on the EE stabilization in 3D space. Accordingly, the experimental results of the NBM and RJM sets are given in Figs. 8 and 9, respectively. Note that the EE motion is defined w.r.t. $\Sigma_{B,t_{0}}$ and the absolute EE position is defined w.r.t. the origin of the optical tracking system. Moreover, we can see clearly that the z-axis and x-axis are coupled in the U-D motion and the x-axis and y-axis are coupled in the F-B motion. The planned joint velocity trajectories (Figs. 8b and 9b) fall in the joint velocity limit $(-0.1,0.1)$ (rad/s). The EE motion (Figs. 8d and 9d) has the opposite trend w.r.t. the human motion (Figs. 8a and 9a). Figs. 8e and 9e show that EE goes back to its initial orientation after the compensation. Figs. 8f and 9f show the compensation performance in Cartesian space.

Two indices are used to evaluate the 1D motion compensation performance.

	$$\displaystyle\bar{\rho}_{E\chi}=\frac{\sum_{i=k}^{n+k}{\rho}_{E\chi}\left(i\right)}{n}$$		(31)
	$$\displaystyle s_{\chi}=\sqrt{\frac{\sum_{i=k}^{n+k}\left({\rho}_{E\chi}\left(i\right)-\bar{\rho}_{E\chi}\right)^{2}}{n-1}}$$		(32)

where ${\rho}_{E\chi}$ represents the absolute EE position obtained from the optical tracking system, $\bar{\rho}_{E\chi}$ represents the mean value of ${\rho}_{E\chi}$, $\chi$ represents x, y, and z coordinates, and $s_{\chi}$ represents the standard deviation of ${\rho}_{E\chi}(i)$. Meanwhile, we defined ${e}_{\chi}=\left|\bar{\rho}_{E\chi}-\bar{\rho}_{EI\chi}\right|$ to represent the mean EE position error, where $\bar{\rho}_{EI\chi}$ is the initial mean value of ${\rho}_{E\chi}$ $(t\leq 3s)$.

In the NBM set, the initial positions are $\bar{\rho}_{EIx}=0.1280$, $\bar{\rho}_{EIy}=-0.5606$, and $\bar{\rho}_{EIz}=1.9089$. For the U-D motion, $\bar{\rho}_{Ez}=1.9377$, $s_{z}=0.0390$, and ${e}_{z}=0.0288$. For the L-R motion, $\bar{\rho}_{Ey}=-0.5291$, $s_{y}=0.0411$, and ${e}_{y}=0.0315$. For the F-B motion, $\bar{\rho}_{Ex}=0.0884$, $s_{x}=0.0191$, and ${e}_{x}=0.0396$. By comparison, in the RJM set, the initial positions are $\bar{\rho}_{EIx}=-0.0188$, $\bar{\rho}_{EIy}=-0.4993$, and $\bar{\rho}_{EIz}=1.9649$. For the U-D motion, $\bar{\rho}_{Ez}=1.9489$, $s_{z}=0.0497$, and ${e}_{z}=0.0160$. For the L-R motion, $\bar{\rho}_{Ey}=-0.4589$, $s_{y}=0.0347$, and ${e}_{y}=0.0404$. For the F-B motion, $\bar{\rho}_{Ex}=-0.0024$, $s_{x}=0.0248$, and ${e}_{x}=0.0164$. Furthermore, we propose a distance-related index to evaluate the 3D motion compensation in Cartesian performance.

$$\displaystyle D_{E}=\frac{\sqrt{e_{x}^{2}+e_{y}^{2}+e_{z}^{2}}}{\sqrt{{e}_{x,max}^{2}+{e}_{y,max}^{2}+{e}_{z,max}^{2}}}$$

(33)

where ${e}_{\chi,max}=max(|\bar{\rho}_{E\chi}-{\rho}_{E\chi}\left(i\right)|)$. In the 3D random motion, $D_{E}$ of the NBE set is $0.3939$ and $D_{E}$ of the RJM set is $0.3709$. Meanwhile, the evaluation results ($\bar{\rho}_{E\chi}$, $s_{\chi}$, and $D_{E}$) of two experiments are illustrated in Fig. 10.

From the experimental results, we can see that the proposed control framework can effectively compensate for the human translation motion. However, the compensation motion sometimes falls behind the user motion due to the joint clamping (for safety) and the sampling frequency of the IMU-based sensory interface (60 Hz). For example, if the previous step compensation hasn’t been finished, the motion compensation in Fig.6 (d-g) still goes backward. From the evaluation results, we can see that both methods have the same order of magnitude in ${e}_{\chi}$ and $D_{E}$ (RJM is slightly better than NBM yet simpler).