Bond Graph Modeling and Kalman Filter Observer Design for an Industrial Back-Support Exoskeleton

Abstract

This paper presents a versatile approach to the synthesis and design of a bond graph model and a Kalman filter observer for an industrial back-support exoskeleton. Actually, the main purpose of developing a bond graph model is to investigate and understand better the system dynamics. On the other hand, the design of the Kalman observer always should be based on a model providing an adequate description of the system dynamics; however, when back-support exoskeletons are considered, the synthesis of a state observer becomes very challenging, since only nonlinear models may be adopted to reproduce the system dynamic response with adequate accuracy. The dynamic modeling of the exoskeleton robotic platform, used in this work, comprises an electrical brushless DC motor, gearbox transmission, torque sensor and human trunk (biomechanical model). On this basis, a block diagram model of the dynamic system is presented and an experimental test has been carried out for identifying the system parameters accordingly. Both the block diagram and bond graph dynamic models are simulated via MATLAB and 20-sim software (bond graph simulation software) respectively. Furthermore, the possibility of employing the Kalman filter observer together with a suitable linear model is investigated. Subsequently, the performance of the proposed Kalman observer is evaluated in a lifting task scenario with the use of a linear quadratic regulator (LQR) controller with double integral action. Finally, the most important simulation results are presented and discussed.

1. Introduction

Considerable research effort has recently been devoted to the development of adequate models and control approaches for wearable robots such as exoskeletons [1]. Generally, an exoskeleton is worn by an operator to provide support and protection from physical activities by generating proper assistive force(s) and torque(s) [2]. Additionally, exoskeletons might be profitably used in such non-industrial areas such as military, medicine and rehabilitation [3]. For military applications, the exoskeletons are designed to enhance the travel abilities and loading capabilities of soldiers; for industrial purposes, the exoskeletons are used to augment an operator’s load-handling potential and avoid and even reduce the occurrence of musculoskeletal damages; in rehabilitation, exoskeletons improve the lost functions and the quality of life of patients with serious disabilities, motor cognitive restrictions [4]. An example of a back-support exoskeleton is illustrated in Figure 1.

Due to the direct interface between the operator and the physical structure there are a number of features in exoskeleton design that are not common, or not so critical, in the design of industrial robots. From the mechanical point of view, motion ranges, safety, user comfort, low inertia and adaptability should be taken into account [5]. On the other hand, controllability, responsiveness, flexible and smooth motion generation should be considered in exoskeleton control schemes [6].

With the intention of applying control schemes to exoskeletons and improving the system performance, a comprehensive dynamic model to represent the behaviour of this robotic platform is needed [7]. Simple dynamic models for exoskeleton robotic platforms such as the biomechanical model in 2D dimensional space [8], the Lagrangian dynamics model, specialized with a serial rigid links chain [9] and the mathematical model based on user’s movements [10] have been presented so far. Basically, precise dynamic model of the robotic platform represents the system behavior more accurately, which results in more appropriate and desirable control performance [11]. In addition, while achieving model accuracy for complex systems, the order of the system increases which results in complicated computations [12]. Modeling of complex systems such as exoskeletons through bond graphs requires fewer mathematical computations and is more adjustable and adaptable with respect to classical conventional dynamic modeling approaches. In this modeling approach, the dynamic systems are exhibited graphically, subsystems interconnected through bonds based on the flow of power in the system. This makes bond graph modeling an appropriate approach for modeling multidisciplinary dynamic systems [13].

On the other hand, one of the main issues for developing exoskeleton systems is the lack of appropriate measuring devices and methods for the direct measurement of all the system states. This issue can limit the use of some control techniques such as model-based controls [14] in terms of the state feedback [15] since model-based controls assume that the system dynamic is known [16]. More recently, Yang, et al. in [17] proposed a novel compliant upper-body exoskeleton platform to reduce spine loading across multiple vertebral joints. Although their robotic control system demonstrates a good robust behaviour, the tracking error is still high due to lack of information from the interaction forces and torques between the human body and exoskeleton. In the literature, several methods have been proposed to measure and estimate the human-robot interaction force/torque in exoskeleton applications. Choi in [18] developed and proposed a compact force sensor system for an exoskeleton platform to measure interaction force/torque. However, using force/torque sensors are generally expensive and cause an increase in the device weight and complexity. Some proposed methods rely on commanded actuation torques (motor currents) to estimate the interaction force/torque such as the work presented in [19]. Distinguishing the damping and elasticity effects of system components for identifying system dynamics is the main problem of these methods. Furthermore, employing EMG sensors to detect the biological signals from muscle activities to estimate interaction force/torque is another popular approach among researchers [20,21]. Nevertheless, EMG sensors are highly sensitive to to changes in skin conductivity skin conductivity, mounting and positioning, muscle fatigue and interactions between muscles. Using inverse dynamics [22], requiring high computational cost, is another method to estimate interaction force/torque. On this basis, the synthesis of a Kalman filter observer for a back-support exoskeleton for estimating interaction force/torque has been proposed and developed as one of the aims of this paper.

In this paper, a dynamic model, employing the bond graph modeling method, rather than the more conventional block diagram approach has been developed, presented and used to better analyze the behavior of the system. In other words, one can easily find out and analyze the effects of elements and variables on the dynamics of an industrial back-support exoskeleton platform. The bond graph modeling technique brings benefits in comparison with the other classical modeling techniques [23]:

• Providing a visual representation of main properties of the system by separating the system formation from dynamic modeling equations;
• Monitoring the coherence and consistency of the topological structure of the dynamic system;
• Providing a hierarchical dynamic model of the system.

Additionally, the synthesis of a Kalman filter observer for a back-support exoskeleton, besides bond graph modeling, is discussed and developed in this work. Kalman filters are usually very sensitive to the modeling errors. However, they can work well in closed-loop system schemes which provide good stability and performance despite large errors in the estimated values. The proposed Kalman observer needs only the measurement of the angular velocity of the electrical motor and the instant value of the torque signal from the output of the transmission section. Hence, the use of other sensors such as strain bridges or accelerometers are not required. The performance of the Kalman state observer has been examined in the simulation an experimental tests of a lifting task scenario by using a linear quadratic regulator (LQR) controller [24] with double integral action. The work, proposed and discussed in this paper, features a high level of originality, rather than the more conventional block diagram approach has been developed, presented and used to better analyze the behavior of the system.

The present paper aims to develop a bond graph model for industrial back-support exoskeletons as an alliterative to classical dynamic modeling methods to describe dynamic behaviors of the system. Furthermore, a general approach for designing a functional Kalman filter observer by means of a linear dynamic model for industrial back-support exoskeletons has been proposed. This formation overcomes the difficulties caused by geometric and technical structural constraints in of back-support exoskeleton dynamic models.

This paper is organized as follows: first, in Section 2, the dynamic model of a back-support exoskeleton is described in terms of block diagram and bond graph modeling. In Section 3, the synthesis of the Kalman State observer is given. The system identification procedure based on the experimental tests is discussed in Section 4. The performance analyses of both the bond graph model and Kalman state observer is discussed in Section 5. Finally, in Section 6, conclusions are drawn.

2. Dynamic Model

In this section, the dynamic model of an industrial back-support exoskeleton with respect to the block diagram and bond graph modeling approach is explained. Theoretically, wearing an exoskeleton results in interaction forces between the user’s body and the robotic structure. These interaction forces can be formulated and modeled through the linear mass-spring-damper system [25]. In this modeling approach, the masses reproduce the inertia properties of different body segments which includes hard and soft tissues. The biomechanical properties of the other segments comprising bones, muscles, tendons, and ligaments are represented by springs and dampers.

2.1. Block Diagram Model

To formulate the the movement of the industrial back-support exoskeleton worn by the operator, a corresponding linear biomechanical model is shown in Figure 2. For the sake of simplicity, a graphical representation of the dynamic model of an industrial back-support exoskeleton is illustrated in Figure 3. The dynamic model comprises a brushless DC motor coupled through the rotor shaft to a gearbox transmission, a torque sensor and a the biomechanical model of the human trunk. Consequently, the dynamic model of the system can be formulated mathematically as:

$$\begin{array}{c}\hfill {\tau }_m-b_f{\dot{\theta }}_m-(J_m+J_g){\ddot{\theta }}_m-\frac{k_g}{r_g}\left((\frac{{\theta }_m}{r_g}-\theta s_{int})\right)-\frac{d_g}{r_g}\left((\frac{{\dot{\theta }}_m}{r_g}-\dot{\theta }{s}_{int})\right)-\frac{b_g}{r_g}({\dot{\theta }}_m)=0\\ \hfill k_g\left((\frac{{\theta }_m}{r_g}-\theta s_{int})\right)+d_g\left((\frac{{\dot{\theta }}_m}{r_g}-\dot{\theta }{s}_{int})\right)-J{s}_{int}\ddot{\theta }{s}_{int}-k_s(\theta s_{int}-\theta s_{ext})-d_s(\dot{\theta }{s}_{int}-\dot{\theta }{s}_{ext})=0\\ \hfill k_s(\theta s_{int}-\theta s_{ext})+d_s(\dot{\theta }{s}_{int}-\dot{\theta }{s}_{ext})-J{s}_{ext}\ddot{\theta }{s}_{ext}-k_{eh}(\theta s_{ext}-{\theta }_h)-d_{eh}(\dot{\theta }{s}_{ext}-{\dot{\theta }}_h)=0\\ \hfill k_{eh}(\theta s_{ext}-{\theta }_h)+d_{eh}(\dot{\theta }{s}_{ext}-{\dot{\theta }}_h)-J_h{\ddot{\theta }}_h-b_h{\dot{\theta }}_h=0\end{array}$$

(1)

The parameters and variables of Equation (1) are defined in

Table 1. Definitions of parameters and variables used in the dynamic model.

Symbol	Quantity	Unit
${\tau }_m$	Torque generated by brush-less DC motor	$[\mathrm{N} \mathrm{m}]$
$b_f$	Friction between brushless DC motor and gearbox	$[\mathrm{N} \mathrm{m} \mathrm{s} {\mathrm{rad}}^{-1}]$
${\theta }_m$	Angular position of brushless DC motor	$\left[\mathrm{rad}\right]$
${\dot{\theta }}_m$	Angular velocity of brushless DC motor	$[\mathrm{rad} {\mathrm{s}}^{-1}]$
${\ddot{\theta }}_m$	Angular acceleration of brushless DC motor	$[\mathrm{rad} {\mathrm{s}}^{-2}]$
$J_m$	Inertia of brushless DC motor	$[\mathrm{k}\mathrm{g} {\mathrm{m}}^2]$
$J_g$	Inertia of gearbox transmission	$[\mathrm{k}\mathrm{g} {\mathrm{m}}^2]$
$k_g$	Gearbox stiffness	$[\mathrm{N} \mathrm{m} {\mathrm{rad}}^{-1}]$
$r_g$	Gearbox ratio	$[-]$
$d_g$	Gearbox damping	$[\mathrm{N} \mathrm{m} \mathrm{s} {\mathrm{rad}}^{-1}]$
$b_g$	Gearbox friction	$[\mathrm{N} \mathrm{m} \mathrm{s} {\mathrm{rad}}^{-1}]$
$\theta s_{int}$	Angular position of torque sensor internal ring	$\left[\mathrm{rad}\right]$
$\dot{\theta }{s}_{int}$	Angular velocity of torque sensor internal ring	$[\mathrm{rad} {\mathrm{s}}^{-1}]$
$\ddot{\theta }{s}_{int}$	Angular acceleration of torque sensor internal ring	$[\mathrm{rad} {\mathrm{s}}^{-2}]$
$J{s}_{int}$	Inertia of torque sensor internal ring	$[\mathrm{k}\mathrm{g} {\mathrm{m}}^2]$
$d_s$	Torque sensor damping	$[\mathrm{N} \mathrm{m} \mathrm{s} {\mathrm{rad}}^{-1}]$
$k_s$	Torque sensor stiffness	$[\mathrm{N} \mathrm{m} {\mathrm{rad}}^{-1}]$
$\theta s_{ext}$	Angular position of torque sensor external ring	$\left[\mathrm{rad}\right]$
$\dot{\theta }{s}_{ext}$	Angular velocity of torque sensor external ring	$[\mathrm{rad} {\mathrm{s}}^{-1}]$
$\ddot{\theta }{s}_{ext}$	Angular acceleration of torque sensor external ring	$[\mathrm{rad} {\mathrm{s}}^{-2}]$
$J{s}_{ext}$	Inertia of torque sensor external ring	$[\mathrm{k}\mathrm{g} {\mathrm{m}}^2]$
$d_{eh}$	Equivalent damping between trunk and exoskeleton	$[\mathrm{N} \mathrm{m} \mathrm{s} {\mathrm{rad}}^{-1}]$
$k_{eh}$	Equivalent stiffness between trunk and exoskeleton	$[\mathrm{N} \mathrm{m} {\mathrm{rad}}^{-1}]$
${\theta }_h$	Angular position of human trunk	$\left[\mathrm{rad}\right]$
${\dot{\theta }}_h$	Angular velocity of human trunk	$[\mathrm{rad} {\mathrm{s}}^{-1}]$
${\ddot{\theta }}_h$	Angular acceleration of human trunk	$[\mathrm{rad} {\mathrm{s}}^{-2}]$
$J_h$	Human trunk inertia	$[\mathrm{k}\mathrm{g} {\mathrm{m}}^2]$
$b_h$	Human hip damping	$[\mathrm{N} \mathrm{m} \mathrm{s} {\mathrm{rad}}^{-1}]$

. By considering the human trunk velocity (${\ddot{\theta }}_h$$[\mathrm{rad} {\mathrm{s}}^{-1}]$ ) and motor torque (${\tau }_m$) as the system inputs, the block diagram corresponding to the dynamic model of the industrial exoskeleton system can be depicted in Figure 4. This block diagram aims to reflect the constitution of the mathematical model. The codes for modeling the dynamic system, developed in the MATLAB environment, are given in Supplementary Materials.

2.2. Bond Graph Model

The bond graph modeling approach is a graphical technique to model and analyze the dynamics of physical dynamic systems from different domains to in a united style. In bond graph modeling approach, the ports of components are connected by bonds that specify the transfer of energy between dynamic system components. These components are described through their energetic behavior. In fact, they can either store or dissipate, supply or absorb, and transform energy reversibly or irreversibly [26]. More information and details about bond graph methodology are given in [27].

Bond graphs causality is a notation format that indicates which side of a bond determines the instantaneous effort and which determines the instantaneous flow. From the graphical point of view, the causality is shown by putting a vertical line close to the element which dominates the flow as illustrated in Figure 5. In dynamic equation formulations that describe the system, causality determines the dependency of variables which correspond to each dynamic modeling element. Consequently, a bond graph can present from different fields and domains in a unified way.

Table 2. Bond graph variable analogies.

Energy Domain	Effort Variable	Flow Variable
Electrical	Voltage	Current
Mechanical	Force	Velocity
Hydraulic	Pressure	Volume flow
Thermal	Temperature	Entropy flow

reports the effort and flow variables of different energy domains.

Based on the block diagram modeling descriptions in the earlier section, the bond graph of the exoskeleton framework can be shown in Figure 6. As can be seen from Figure 2, the exoskeleton robotic platform includes four main parts: brushless DC motor, gearbox transmission, torque sensor and human trunk. Consequently, $S_e$ corresponds to the output torque (${\tau }_m$) of the brushless DC motor in Figure 6. The output torque of the brushless DC motor (${\tau }_m$) drives the motor’s shaft with an inertia of $J_m$ as well as the gearbox transmission with input angular inertia of $J_g$ and it also overcomes the mechanical friction $R_1$:${\tau }_f$. Because these elements are driven with the same velocity, the first (from the left to right) is a 1-type junction. In this work, the shaft of the brushless DC motor is assumed to be stiff.

The gearbox transmission is modeled as transformer ($TF$) and its ratio number is represented as the transformer modulus ($r_g$). The transformed velocity by the gearbox transmission which is effected by the gearbox friction ($b_g$) is given by the 1-type junction right to the $TF$ element.

The speed difference between the gearbox transmission and the internal ring of the torque sensor inertia ($I_3$:${J_S}_{int}$) drives the gearbox stiffness ($C_1$:$\frac{1}{k_g}$) and gearbox damping parameter ($R_3$:$d_g$) through a 0-type junction. The 0-type junctions generate the difference of speed maintaining a constant effort, while the second and third 1-type junction (from left to right) feed the speed differences to both $C_1$:$\frac{1}{k_g}$ and $R_3$:$d_g$ and splits the effort between both elements. In a similar way, the speed difference between the internal and external ring of the torque sensor ($I_4$:${J_S}_{ext}$) drives the torque sensor stiffness ($C_2$:$\frac{1}{k_s}$) and damping ($R_4$:$d_s$) toward a 0-type junction. As much as the torque acts on the external ring of the torque sensor, the corresponding effort is applied to $I_4$:${J_S}_{ext}$ via an 1-type junction.

Regarding to the block diagram model of the system described in Figure 3, the human trunk contains the equivalent damping ($R_5$:$d_{eh}$) and stiffness ($C_3$:$\frac{1}{k_{eh}}$) between the trunk and exoskeleton and also human trunk inertia ($I_5$:$J_h$) and damping ($R_6$:$b_h$). Furthermore, the trunk velocity (${\dot{\theta }}_h$) is considered to be a system input. $R_6$:$b_h$ and $I_5$:$J_h$ have the same velocity (${\dot{\theta }}_h$), while $C_3$:$\frac{1}{k_{eh}}$ and $R_5$:$d_{eh}$ are fed by the difference in the speed between $I_4$:${J_S}_{ext}$ and $I_5$:$J_h$. Therefore, the human trunk inertia and damping are connected to the flow source (trunk velocity) through a 1-junction (5th from left to right).

3. Kalman Filter Observer Design

Basically, a state observer is a system that models a real dynamic system in order to provide an estimation of its internal state variables, given measurements of the input and output of the real dynamic system [28]. In our application, there is no direct measurement for the torque on the human body which is defined as the ratio of the applied torque to the deformation of the trunk muscles. Hence, there is a need for the state observer to obtain an estimation of the torque on the human body. In this section, the design of a Kalman state estimator [29] for a back-support exoskeleton has been explained. The first step to implement the Kalman state observer in our application is to take the state vector into account as:

$$\begin{array}{c}\hfill x={\left[{\theta }_m{\dot{\theta }}_m\theta s_{int}\dot{\theta }{s}_{int}\theta s_{ext}\dot{\theta }{s}_{ext}{\theta }_h\right]}^T\end{array}$$

(2)

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\dot{\theta }}_m\\ {\ddot{\theta }}_m\\ \dot{\theta }{s}_{int}\\ \ddot{\theta }{s}_{int}\\ \dot{\theta }{s}_{ext}\\ \ddot{\theta }{s}_{ext}\\ {\dot{\theta }}_h\end{array}\right]=\left[\begin{array}{ccccccc}0& 1& 0& 0& 0& 0& 0\\ \frac{-k_g}{(J_m+J_g)r_g^2}& \frac{-(b_f.r_g^2+b_g+d_g)}{(J_m+J_g)r_g^2}& \frac{k_g}{(J_m+J_g)r_g}& \frac{d_g}{(J_m+J_g)r_g}& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0\\ \frac{k_g}{J{s}_{int}.r_g}& \frac{d_g}{J{s}_{int}.r_g}& \frac{-(k_g+k_s)}{J{s}_{int}}& \frac{-(d_g+d_s)}{J{s}_{int}}& \frac{k_s}{J{s}_{int}}& \frac{d_s}{J{s}_{int}}& 0\\ 0& 0& 0& 0& 0& 1& 0\\ 0& 0& \frac{k_s}{J{s}_{ext}}& \frac{d_s}{J{s}_{ext}}& \frac{-(k_{eh}+k_s)}{J{s}_{ext}}& \frac{-(d_{eh}+d_s)}{J{s}_{ext}}& \frac{k_{eh}}{J{s}_{ext}}\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right]\times \left[\begin{array}{c}{\theta }_m\\ {\dot{\theta }}_m\\ \theta s_{int}\\ \dot{\theta }{s}_{int}\\ \theta s_{ext}\\ \dot{\theta }{s}_{ext}\\ {\theta }_h\end{array}\right]\\ \hfill +{\left[\begin{array}{c}0\frac{1}{J_m+J_g}00000\end{array}\right]}^T{\tau }_m+{\left[\begin{array}{c}00000\frac{d_{eh}}{J{s}_{ext}}1\end{array}\right]}^T{\dot{\theta }}_h\end{array}$$

(3)

Equation (1) can be reorganized and expressed in the form of a matrix representation (Equation (3)). Based on this, by taking process and measurement noises into the formulation, a more compact form of Equation (3) can be expressed in the form of state-space equations as follows:

$$\left\{\begin{array}{c}{\dot{x}}_1=A_1{x}_1+A_3{x}_2+B_m{\tau }_m+B_{{\dot{\theta }}_{h_1}}{\dot{\theta }}_h+w_{{\tau }_m}\hfill \\ {\dot{x}}_2=B_{{\dot{\theta }}_{h_2}}{\dot{\theta }}_h+w_{{\dot{\theta }}_h}\hfill \\ y_v=C_{v}x+v_v\hfill \\ y_t=C_{t}x+v_t\hfill \end{array}\right.$$

(4)

where ${\dot{x}}_1$ includes the first six states of $\dot{x}$ and ${\dot{x}}_2$ is equal to ${\dot{\theta }}_h$. $A_1\in R^{6\times 6}$ comprising the first six rows and columns of A, $A_3\in R^{6\times 1}$ containing the first six rows of the last column of A , $B_m\in R^{6\times 1}$ consisting the first six rows of $B_m$, $B_{{\dot{\theta }}_{h_1}}\in R^{6\times 1}$ including the first six rows of $B_{{\dot{\theta }}_h}$ and $B_{{\dot{\theta }}_{h_2}}\in R^{1\times 1}$ and is equal to 1.

The measuring output $y_v\in R^{1\times 1}$ and $y_t\in R^{1\times 1}$ consists of the measured signal from the brushless DC motor velocity and the signal from the torque sensor respectively by setting:

$$\begin{array}{c}\hfill C_v=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\end{array}\right]\\ \hfill C_t=\left[\begin{array}{cccccc}0& 0& k_s& d_s& -k_s& -d_s\end{array}\right]\end{array}$$

(5)

A sufficient condition for designing a Kalman observer is is satisfying the observability condition of the system. On this basis and by considering $A^{\prime }=\left[\begin{array}{cc}{A}_1& A_3\\ 0_{6\times 1}& 0\end{array}\right]$ and $C^{\prime }=\left[\begin{array}{c}{C}_v\\ C_b\end{array}\right]$, it can be proved that the pair ($A^{\prime },C^{\prime }$) fulfils the observability requirements for existence of an asymptotically stable Kalman observer. Additionally, the torque on the human body that should be estimated in our system is presented by ${\mathbf{y}}_{\mathbf{b}}$ and can be formulated as a system output:

$$y_b=C_b{x}_1-k_{eh}{x}_2-d_{eh}{\dot{\theta }}_h$$

(6)

in which

$$\begin{array}{c}\hfill C_b=\left[\begin{array}{cccccc}0& 0& 0& 0& k_{eh}& d_{eh}\end{array}\right]\end{array}$$

(7)

the observer model can be expressed by the following system of equations:

$$\left\{\begin{array}{c}\widehat{\dot{x}}=\widehat{A}\widehat{x}+\left[B_m{B}_{{\dot{\theta }}_h}\right]\left[\begin{array}{c}{\tau }_m\\ {\dot{\theta }}_h\end{array}\right]u+L(y-\widehat{y})\hfill \\ \widehat{y}=C^{\prime }\widehat{x}\hfill \end{array}\right.$$

(8)

where $\widehat{A}=A^{\prime }-L{C}^{\prime }$ and is asymptotically stable. L is Kalman gain and is chosen to minimize the mean square error between the estimated and the actual state variables. It can be proven that $L=P{C}^{\prime T}{R}^{-1}$, where P is symmetric and a positive semi-definite solution of Riccati equation.

$$A^{\prime }P+P{A}^{\prime T}-P{C}^{\prime T}{R}^{-1}{C}^{\prime }P+Q=0$$

(9)

where $Q\in R^{7\times 7}$ and $R\in R^{2\times 2}$ are the process and measurement noise covariance matrices. Additionally, $\widehat{x}$ represents estimated states while the estimated output is denoted by $\widehat{y}$ in Equation (7). The flowchart in Figure 7 shows the overall algorithm for the Kalman filter developed in this work.

4. System Identification

Friction and damping parameters have a fundamental role in the system dynamic and characterizing these parameters improves the model accuracy.. The value of $b_f$ was estimated from the experimental tests from the work in [7] using the same actuator and it is given by $b_f=0.48\times {10}^{-3}$$[\mathrm{N} \mathrm{s} {\mathrm{rad}}^{-1}]$.

Some system parameters can be identified from components data-sheets provided by manufacturers; however, the information related to some parameters such as $b_g$ and $d_g$ are not reported anywhere. Hence, the experimental setup used for identification of $b_g$ and $d_g$ (shown in Figure 8) developed. This experimental setup comprises the same components used in our exoskeleton. In this setup, an inertia of $0.0491$$[\mathrm{k}\mathrm{g} {\mathrm{m}}^2]$ is attached to the system in place of the human inertia ($J_h$). Therefore in the model used for the identification test $J_h$$=0.0491$$[\mathrm{k}\mathrm{g} {\mathrm{m}}^2]$. The information of main components is reported in

Table 3. Main component used in the experimental setup.

Component	Model
Brushless DC Motor	Maxon-EC 90 flat, 90 Watt
Gearbox	Harmonic Drive, SHD-25-160-2SH-SP
Torque Sensor	ME, TS110a ± 100 Nm

.

To identify the aforementioned unknown system parameters, we compare an experimental open-loop frequency sweep with the dynamic model presented in Section 2. To implement the experimental test, a current chirp reference signal with an amplitude of $2.3$ $\left[\mathrm{A}\right]$, which is equivalent to a torque with the amplitude of $0.25$ $\left[\mathrm{N}\right]$ and the starting frequency ($f_0$) of $0.1$ $\left[\mathrm{Hz}\right]$ and the final frequency ($f_f$) of 130 $\left[\mathrm{Hz}\right]$, was applied to the brushless DC motor. The bode plots extracted from the collected experimental data are shown in Figure 9. Correspondingly, the input data is the torque produced by the motor and the output is the torque signal measured by the torque sensor. By setting $b_g=0.65$ $[\mathrm{N} \mathrm{s} {\mathrm{rad}}^{-1}]$ and $d_g=6$ $[\mathrm{N} \mathrm{s} {\mathrm{rad}}^{-1}]$ a very high accuracy mathematical model corresponding to the experimental data can be obtained as is shown in Figure 9.

To obtain a more reliable dynamic model from the system, noise effects have been taken into consideration on system inputs (process) and outputs (measurements). In this work, $w_{{\tau }_m}$$\left[\mathrm{N}\right]$ and $w_{{\dot{\theta }}_h}$$[\mathrm{rad} {\mathrm{s}}^{-1}]$ represent noise parameters on the motor torque input (${\tau }_m$) and human velocity input (${\dot{\theta }}_h$) respectively. On the other hand, noise parameters on motor encoder and torque sensor signals are given by $v_v$ $[\mathrm{rad} {\mathrm{s}}^{-1}]$ and $v_t$ $\left[\mathrm{N}\right]$ respectively. The range used for these noise parameters are chosen based on the robot components as: $|w_{{\tau }_m}|\le 0.1\times {10}^{-3}$ $\left[\mathrm{N}\right]$, $|w_{{\dot{\theta }}_h}|\le 0.15\times {10}^{-3}$ $[\mathrm{rad} {\mathrm{s}}^{-1}]$, $|v_v|\le 0.1\times {10}^{-6}$ $[\mathrm{rad} {\mathrm{s}}^{-1}]$ and $|v_t|\le 35\times {10}^{-3}$ $\left[\mathrm{N}\right]$.

In addition, it should be mentioned that the values for $k_{eh}$ and $b_h$ are extracted and estimated for a person with a weight of 75 $\left[\mathrm{k}\mathrm{g}\right]$ and a height of $1.75$ $\left[\mathrm{m}\right]$ as 1125 $[\mathrm{N} {\mathrm{rad}}^{-1}]$ and $0.75$$[\mathrm{N} \mathrm{s} {\mathrm{rad}}^{-1}]$ respectively based on the work presented in [30]. Furthermore, the friction torque between the motor and gearbox transmission ($b_f$), the damping coefficient between the human trunk and robotic platform (exoskeleton) ($d_{eh}$) and torque sensor damping ($d_s$) are neglected as their values are so small with respect to those associated with the motor/gearbox and they can easily be compensated.

5. Simulation Results

5.1. Comparison between Block Diagram and Bond Graph Models

The simulation results from the block diagram and bond graph models are presented and discussed in this section. While the block diagram dynamic model of the exoskeleton robotic system is simulated via MATLAB software, the simulation of the bond graph modeling is carried out via 20-sim software [7]. The results from both simulation softwares are plotted with MATLAB.

The exoskeleton system, presented and discussed in previous sections, has been fed with an input torque of $0.1$ $\left[\mathrm{N}\right]$ and a human trunk velocity of $0.05$ $[\mathrm{rad} {\mathrm{s}}^{-1}]$ with a delay of $0.1$ $\left[\mathrm{s}\right]$. The simulations are carried out based on the torque sensed by the torque sensor as an effort variable in the terminology of the bond graph modeling and also the rotational velocity of the brushless DC motor shaft (gearbox transmission input velocity) is a flow variable in the terminology of the bond graph modeling. However, the simulation tests can be extended and implemented for all the system variables.

A comparison between the torque applied on the torque sensor through results from block diagram and bond graph models is demonstrated in Figure 10. From Figure 10, it can be inferred that the results from the bond graph model for the torque implemented on the torque sensor is very similar and close to that from the block diagram; however, the accuracy is slightly reduced as the robotic platform (human trunk) moves. The whole average error for the simulation test is around $0.09$ $\left[\mathrm{N}\right]$ over 1 $\left[\mathrm{s}\right]$.

In addition, the difference between the rotational velocity of the brushless DC motor shaft through results from block diagram and bond graph models is shown in Figure 11. It can be concluded from Figure 11 that the two models are very close to each other and the overall average error is about $0.07$ $[\mathrm{rad} {\mathrm{s}}^{-1}]$ during the first second. Accordingly, the error at the start of the movement slightly increases, but then it starts to approach zero over time.

Finally, an excellent agreement between the block diagram model developed in MATLAB and bond graph model developed in 20-sim can be inferred from both Figure 10 and Figure 11. It should be pointed out that using different numerical differential equation solvers with diverse configurations by MATALB and 20-sim is the main reason of the small errors between the simulation results [7].

5.2. Implementation of Kalman Filter Observer and Performance Analysis

The theory described in the forgoing sections means that both the dynamic model of the exoskeleton and Kalman observer can be implemented. The performance of the observer has been evaluated by comparing the values of the state variables computed by the simulator with the observer estimates on a test case using of a linear quadratic regulator (LQR) controller [31] with integral action. It is important to note that the synthesis and application of the Linear-Quadratic-Integral control is not described in detail as it is not pertinent to the Kalman observer design and performance assessment. All the simulation results have been executed by means of MATLAB.

The following state-space equations for the LQR controller with integrator action are taken into consideration:

$$\underset{{\dot{x}}_i^{\prime }}{\underbrace{\left[\begin{array}{c}\dot{x}\\ \dot{z}\end{array}\right]}}=\underset{A_i^{\prime }}{\underbrace{\left[\begin{array}{cc}{A}_1& 0_{6\times 1}\\ -C_b& 0_{1\times 1}\end{array}\right]}}\underset{x_i^{\prime }}{\underbrace{\left[\begin{array}{c}x\\ z\end{array}\right]}}+\underset{B_i^{\prime }}{\underbrace{\left[\begin{array}{c}{B}_m\\ 0\end{array}\right]}}{\tau }_m+\underset{B_h^{\prime }}{\underbrace{\left[\begin{array}{c}{B}_{{\dot{\theta }}_h}\\ 0\end{array}\right]}}{\dot{\theta }}_h+\left[\begin{array}{c}0\\ 1\end{array}\right]{\tau }_{ref}$$

(10)

where ${\tau }_{ref}$ is the reference torque of the closed-loop control system and z is the auxiliary state corresponding to $y_b$. Thereby, the optimal state-feedback control law is obtained in the form:

$$u^{\prime }=-F_i^{\prime }{x}_i^{\prime }=-F_{ix}^{\prime }x-F_{iz}^{\prime }z$$

(11)

where $F_i^{\prime }$ represents the optimal feedback matrix gain, $F_{ix}^{\prime }$ is the proportional matrix gain and $F_{iz}^{\prime }$ indicates the integral gain. ${\mathbf{Q}}_{\mathbf{c}}^{\prime }\in {\mathbb{R}}^{7\times 7}$ and ${\mathbf{R}}_{\mathbf{c}}^{\prime }\in {\mathbb{R}}^{1\times 1}$ are positive definite symmetric weighing matrices of the LQR cost function corresponding to the states from the Kalman observer and the human trunk velocity. It is worth mentioning that the states from the observer contains estimated values corresponding to ${\theta }_m,{\dot{\theta }}_m,\theta s_{int},\dot{\theta }{s}_{int},\theta s_{ext},\dot{\theta }{s}_{ext}$ and the estimated torque on the human body.

Varying $Q_c^{\prime }$ and $R_c^{\prime }$ matrices cause changes in the location of the closed-loop poles which effects the performance of the system. While the large value of $R_c^{\prime }$ causes smaller control inputs and thereby larger values of the states, a faster convergence of the states could be obtained by larger elements of $Q_c^{\prime }$. The following consideration are used to achieve a proper feedback matrix:

$$\begin{array}{cccccc}\hfill |{\theta }_m|& \le 0.001\left[\mathrm{rad}\right],\hfill & \hfill |{\dot{\theta }}_m|& \le 0.002[\mathrm{rad} {\mathrm{s}}^{-1}],\hfill & \hfill |\theta s_{int}|& \le 0.0001\left[\mathrm{rad}\right]\hfill \\ \hfill |\dot{\theta }{s}_{int}|& \le 0.0002[\mathrm{rad} {\mathrm{s}}^{-1}],\hfill & \hfill |\theta s_{ext}|& \le 0.0001\left[\mathrm{rad}\right],\hfill & \hfill |\dot{\theta }{s}_{ext}|& \le 0.0002[\mathrm{rad} {\mathrm{s}}^{-1}]\hfill \\ \hfill |e_{\tau }|& \le 15\times {10}^{-6}[\mathrm{N} \mathrm{m}],\hfill & \hfill |{\tau }_m|& \le 0.56[\mathrm{N} \mathrm{m}]\hfill & \hfill & \end{array}$$

where $e_{\tau }$ is equal to the difference between estimated and reference torque on the human body. To reject the ramp trajectories arising from the human trunk position during the constant human velocity, the LQR requires a double integrator. Hence, the following formulation is derived:

$$\begin{array}{c}\hfill A_i=\left[\begin{array}{cc}{A}_1& 0_{6\times 2}\\ -C_b& 0_{1\times 2}\\ 0000& 0010\end{array}\right]\mathrm{and}{B}_i=\left[\begin{array}{c}{B}_m\\ 0\\ 0\end{array}\right]\end{array}$$

(12)

The LQR optimal control law of Equation (12) can be written as:

$$u=-F_{ix}^{\prime }x-F_{iz}^{\prime }z-F_{iz}{\int }_0^{t}zd\tau$$

(13)

in which t represents the instantaneous time and $\tau$ is the variable of the integration associated with the time.

The controllability condition of the pair ($A_i,B_i$) is an adequate condition for the presence of a state feedback gain $F_i$ in a way that ($A_i-B_i{F}_i$) be asymptotically stable. $F_i$, the controller matrix, can be calculated by minimizing the quadratic performance index of the LQR controller for Equation (12). Figure 12 demonstrates the block diagram of the closed-loop control system. In Figure 12, $F_{iz}$ is the gain related to the integrator term corresponding to the torque on the human trunk (last array of $F_i$) while the estimated torque on the human trunk is indicated by $\widehat{y_b}$.

A simulation to evaluate the observer performance has been carried out. The main assumption of this simulation is that the exoskeleton works in assistive mode. That is, the human trunk moves with a constant velocity (${\dot{\theta }}_h=0.1$ $[\mathrm{rad} {\mathrm{s}}^{-1}]$) and an assistive torque is applied on the human body by the exoskeleton. This simulated scenario is archived by applying a step reference as the system input with an amplitude of ${\tau }_{ref}=$ 10 $[\mathrm{N} \mathrm{m}]$ and a delay of $0.1$ $\left[\mathrm{s}\right]$ for a time period of $0.5$ $\left[\mathrm{s}\right]$.

It is worth noting that the resulting $Q_c^{\prime }$ and $R_c^{\prime }$ matrices and weighting matrices (${\mathbf{Q}}_{\mathbf{c}}\in {\mathbb{R}}^{8\times 8}$ and ${\mathbf{R}}_{\mathbf{c}}\in {\mathbb{R}}^{1\times 1}$) associated with LQR related to Equation (12) are:

$$\begin{array}{c}\hfill Q_c^{\prime }=diag({10}^5,{10}^5,{10}^7,{10}^7,{10}^7,{10}^7,5\times {10}^{13})\mathrm{and}{R}_c^{\prime }={10}^{12}\end{array}$$

(14)

$$\begin{array}{c}\hfill Q_c=diag({10}^5,{10}^5,{10}^7,{10}^7,{10}^7,{10}^7,5\times {10}^{13},5\times {10}^{15})\mathrm{and}{R}_c={10}^{12}\end{array}$$

(15)

The final state feed-backs are extracted as follows:

$$\begin{array}{c}\hfill F_{ix}^{\prime }=\left[1.000.030.6400.020\right]\end{array}$$

(16)

$$\begin{array}{c}\hfill F_{iz}^{\prime }=[-3.16]\mathrm{and}{F}_{iz}=[-70.71]\end{array}$$

(17)

Essentially, controllers based on state feedback which are implemented through state observers can have disappointing properties regarding relative stability even when the state feedback and state observers demonstrate good and robust stability characteristics [32]. On this basis, stability margins have been evaluated by Nyquist plot (Figure 13 and Figure 14) of the loop gain, from the motor plant input (point $P_{in}$ in Figure 12) to the computed torque signal from the controller (point $P_{out}$ in Figure 12). We can see that the plot in Figure 13, does not encircle the critical point ($-1+0j$), so the system is stable. Figure 14 shows the Nyquist plot in the neighborhood of the critical point. The nominal design has a minimum gain margin of approximately $8.46$$\left[\mathrm{dB}\right]$ and minimum phase margin of about 59 $\left[{}^{\circ }\right]$ with a delay margin of around $0.007$$\left[\mathrm{s}\right]$ which indicates that the minimum sampling frequency of the control unit should be almost 143 $\left[\mathrm{Hz}\right]$ to guarantee stability.

The bode plot of the transfer function from the input disturbance (human trunk velocity) to the estimated torque on the human body is shown in Figure 15 to assess the disturbance rejection bandwidth. This determines how quickly the effects correlated with the human trunk motion decline to an acceptable level, for the exoskeleton platform. By considering the bandwidth as the first frequency in which the gain descends to below $-3 \mathrm{dB}$ of its DC value, it can be inferred that the bandwidth for tracking disturbance (human trunk velocity) rejection is almost 52 $\left[\mathrm{Hz}\right]$; nevertheless, it can be concluded that the human motion with the frequency range of 4∼6.2 $\left[\mathrm{Hz}\right]$ has the most destructive effect on the estimated torque.

The first evidence of the observer’s effectiveness is reported in Figure 16 where a comparison is made between the step responses of the closed-loop control system using the proposed state observer and the linear dynamic model of the system. As it can be inferred from Figure 16, the step response of the observer and the linear model are very similar. Figure 17 shows the error on $\delta y_b=y_b-\widehat{y_b}$ that is, the error between the linear model and the estimator in terms of torque on the human body. This is very small with the mean-squared error of about $0.1\times {10}^{-3}$ $\left[\mathrm{N}\mathrm{m}\right]$ and approaches zero over time. Another interesting finding from Figure 16 is the undershoot at the beginning of the human motion. This occurs because of the phase delay between the human trunk motion velocity and the estimated torque (non-minimum phase system).

The DC motor velocity is shown by Figure 18 and Figure 19. While Figure 18 shows a comparison between the measured and estimated values of the motor velocity, the error between these values ($\delta {\dot{\theta }}_m$) over time is reported in Figure 19. Regarding these figures, the mean-squared error is around $3\times {10}^{-6}$ $[\mathrm{rad} {\mathrm{s}}^{-1}]$ which implies a good agreement between the measured and estimated variables.

Furthermore, according to Figure 20, the overall torque produced by the brushless DC motor is in the defined torque limit range ($\pm 0.56$ $[\mathrm{N} \mathrm{m}]$).

6. Conclusions

In this work, the synthesis and design of a bond graph model and an effective Kalman filter observer for an industrial back-support exoskeleton platform are investigated and presented. Both the bond graph model and Kalman observer, presented in this work, can be extended and implemented on all types of active wearable robots.

After the development the dynamic equations of the robotic platform and extraction of the bond graph model and Kalman state observer of the system, the validation procedure is carried out. The bond graph model, corresponding to an industrial back-support exoskeleton, is validated to allow comparison comparison between the results from bond graph approach and those from the block diagram dynamic model according to the torque sensed by the torque sensor and the velocity of the rotational brushless DC motor shaft. It can be inferred that the responses of the two dynamic models lead to very similar trends, as the errors (the difference between the responses) are quite small. On the other hand, the proposed Kalman filter observer makes use of an analytical dynamic model and leads to the synthesis of an observer with a structure combining the features of an asymptotic Kalman estimator and a linear system. This structure overcomes the difficulties arising from measuring directly the human-robot interaction torque. The results of the observer, achieved from the simulation environment through a Linear-Quadratic-regulator (LQR) control with double integrator action and the full state feedback, illustrate an accurate estimation a very low estimation error. The availability of an effective and computationally efficient state observer represents an essential starting point for the synthesis and experimental validation of model-based control schemes based on state vector feedback.

The main challenge of the presented work is to derive a dynamic model for the human body. There is, of course, a relatively wide variation in the characteristics of the human body and current approaches and data are not yet adequate in these areas.

Future work will consider the experimental validation of both the proposed bond graph model and Kalman filter observer on a real-time hardware by imposing a bound to the frequency content of the simulated system dynamics.