Simulation system, method and computer-readable medium for human augmentation devices

[0001] This application claims priority of U.S. provisional applications No. 60/330,689 filed on Oct. 29, 2001 and No. 60/333,753 filed on Nov. 29, 2001.

[0002] The present invention relates to a simulation system, a simulation method and a computer-readable medium for analysis, and synthesis of human motion under partial assist from powered augmentation devices.

[0003] Effective usage of human assistive systems or augmentation devices for restoration or enhancement of motor function is an important area of research in rehabilitation and performance enhancement. A partial list of important and desired features of an effective assistive system include: (1) a decrease in energy rate and cost with respect to able-bodied subjects performing the same task, (2) minimum disruption and maximum comfort of normal activities when employing the assistive system, and (3) practicality.

[0004] The third requirement considers the ease of wearing such a device and its power consumption needs. These requirements and available technology have led to the development of externally powered orthoses and prostheses that interface directly or indirectly with the human neuromuscular system. Although significant progress has been made in meeting many of the requirements needed for development of practical human assist devices (Popovic, D., Externally Powered and Controlled Orthotics and Prosthetics. The Biomedical Engineering Handbook, Editor Bronzino, J. D., 2nd ed. Vol. 2, Chapter 142, 2000), realization of such systems for daily applications is still in its infancy. The complexity of the central nervous system (CNS) control and the interface between voluntary control and external artificial control are still challenging, unanswered questions.

[0005] At Honda's Wako Research Center, a mechanically powered walking assist prototype system was recently unveiled (Katoh and Hirata, The Concept of a Walking Assistance Suit, Welfare Engineering Symposium, The Japan Society of Mechanical Engineers, August 2001). The target application is to help the elderly and disabled people to either execute daily tasks they could not previously perform, or use less physical exertion than they currently employ for these tasks. The tasks considered include walking, lifting, sitting/standing, and climbing stairs. Two important and challenging questions to consider in the implementation of the Honda prototype and similar human augmentation systems include: 1) analysis and monitoring of biomechanical as well as physiological quantities which cannot be readily measured and 2) the synthesis of an active control which can safely and effectively augment voluntary control. By developing computational methods to study these issues, future performance of human augmentation devices can be studied through simulation, without constraints imposed by hardware implementations of current technology. Simulation studies also enable us to estimate physiological quantities that cannot be easily measured, including muscle forces, joint forces, and energetics of motion. We can simulate effects of aging, predict muscular activity, estimate muscle fatigue and capacity, and detect potential dangerous physiological conditions. It should be mentioned that the exclusive use of simulation is not a substitute for eventual testing on live human subjects. However, an accurate subject-specific simulation allows control algorithms to be designed and refined for the walking assist device. This is especially relevant in our target user population because they already have existing health constraints.

[0006] U.S. Pat. No. 6,152,890 discloses a device and a method for the recording, presentation and automatic classification of biomechanical variables measured on a freely moving test person during a work shift.

[0007] Japanese patent publication unexamined No. 2000-249570 discloses a method for generating human kinematic data.

[0008] Gruber, K. et. al., 1998. A comparative study of impact dynamics: wobbling mass model versus rigid body models. Journal of Biomechanics 31, 439-444 discloses inverse dynamics model used to simulate the human body.

[0009] However, any of the above documents does not deal with analysis and synthesis of human motion under assist from powered augmentation devices.

[0010] Accordingly, what is needed is a system and a method for analysis and synthesis of human motion under assist from powered augmentation devices.

[0011] According to one aspect of the invention, a simulation system is provided for a combined musculoskeletal and augmentation device system including segments and joints connecting the segments. The simulation system comprises a dynamics model of the combined musculoskeletal and augmentation device system and an augmentation device controller for control of the augmentation device. The simulation system further comprises an inverse dynamics model for the musculoskeletal and augmentation device system and a muscle force and muscle capacity module for checking and adjusting the computed torques. The dynamics model of the combined musculoskeletal and augmentation device system receives feasible computed torques at the joints as inputs and delivers simulated kinematic data of the segments as outputs. The augmentation device controller for control of the augmentation device, receives the simulated kinematic data as inputs and delivers assist torques as outputs. The inverse dynamics model for the musculoskeletal and augmentation device system, receives the simulated kinematic data, desired kinematic data of the segments and the assist torques as inputs and delivers the computed muscle torques and net joint torque as outputs. The muscle force and muscle capacity module for checking and adjusting the computed torques, receives the computed muscle torques as inputs and delivers feasible computed torques as outputs after making adjustments to the computed torques.

[0012] According to another aspect of the invention, a method is provided for simulating a combined musculoskeletal and augmentation device system including segments and joints connecting the segments. The method comprises the steps of computing assist torques of the augmentation device, based on simulated kinematic data and computing torques based on the simulated kinematic data, desired kinematic data of the segments and the assist torques. The method further comprises the steps of checking and adjusting the computed torques and computing the simulated kinematic data of the segments based on the computed torques at the joints.

[0013] According to another aspect of the invention, a computer readable medium containing a program for simulating a combined musculoskeletal and augmentation device system including segments and joints connecting the segments, is provided. The program comprises instructions of computing assist torques of the augmentation device, based on simulated kinematic data and of computing torques based on the simulated kinematic data, desired kinematic data of the segments, and the assist torques. The program further comprises instructions of checking and adjusting the computed torques and of computing the simulated kinematic data of the segments based on computed torques at the joints.

[0014] According to an embodiment of the invention, muscle forces are deduced from the computed torques, compared with maximum force limits and adjusted if the muscle forces exceed limits, to obtain feasible torques.

[0015] According to another embodiment of the invention, muscle forces with and without the assist torques are compared in order to asses whether the assist torque control helps or hinders motion and if the assist torque control hinders motion the muscle forces are adjusted and feasible joint torques are computed.

[0016] According to another embodiment of the invention, muscle forces with and without the assist torques are compared in order to asses whether the assist torque control helps or hinders motion and if the assist torque control hinders motion the assist torque control law is adjusted to ensure that feasible joint torques are computed.

[0017] According to another embodiment of the invention, muscle forces are deduced based on a static optimization criterion in which a sum of muscle activation squared is minimized.

[0018] According to another embodiment of the invention, modified accelerations of kinematic data are obtained through non-linear position and velocity feedback from the simulated kinematic data.

[0019] According to another embodiment of the invention, the kinematic data include position data, velocity data and acceleration data and estimates of kinematic data are computed, through non-linear feedback based on desired acceleration data, error between simulated position data and desired position data and error between simulated velocity data and desired velocity data.

[0020] According to another embodiment of the invention, the kinematic data include position data, velocity data and acceleration data and estimates of kinematic data are computed, through non-linear feedback based on error between simulated position data and desired position data and/or error between simulated velocity data and desired velocity data.

[0021] According to another embodiment of the invention, computed reaction forces under the segments contacting the ground are obtained based on the feasible computed torques and the simulated kinematic data.

[0022] According to another embodiment of the invention, gravity compensation control algorithm is employed, in which the assist torques are obtained to reduce the computed muscle force by artificially compensating for the forces due to gravity.

[0023] According to another embodiment of the invention, change in the computed torques, due to compensation for gravity is obtained, using coordinates of the center of the mass of the segments.

[0024] According to another embodiment of the invention, the coordinates of the center of the mass of the segments, are obtained from measurements of joint angles and segment lengths.

[0025] According to another embodiment of the invention, change in the computed torques, due to compensation for gravity, is obtained using measured reaction forces under the feet.

[0026] According to another embodiment of the invention, the feedback gains are selected to produce the fastest possible non-oscillatory response.

[0027] FIG. 1 a biped system having five degrees of freedom in the sagittal plane with intermittent ground contact during double support, single support, and air-born phase;

[0028] FIG. 2 is a system model description with intermittent contact of left and right feet with the ground;

[0029] FIG. 3 is an inverse dynamics controller with position and velocity feedback for calculation of torques that when applied to a system model, will track and reproduce the desired kinematic data;

[0030] FIG. 4 is a muscle force and muscle capacity module;

[0031] FIG. 5 is a block-diagram of the integrated simulation system;

[0032] FIG. 6 is a flowchart illustrating a simulation process according to an aspect of the present invention;

[0033] FIG. 7 is a simulation of the joint angles during a squatting maneuver without employing the desired accelerations (a0), in which nearly perfect tracking of desired kinematic trajectories is illustrated;

[0034] FIG. 8 is a simulation of the joint torques during a squatting maneuver without employing the desired accelerations (a0), in which the proposed method using nonlinear feedback (NLF) produces nearly identical joint torque estimates as compared to the ground truth (ideal) joint torques obtained by a noise free inverse dynamics computation; and

[0035] FIG. 9 is a simulation of the horizontal and vertical ground reaction forces during a squatting maneuver without employing the desired accelerations (a0), in which the proposed method using nonlinear feedback (NLF) produces nearly identical ground reaction estimates as compared to the ground truth (ideal) ground reaction forces obtained by an Iterative Newton Euler inverse dynamics procedure.

[0036] The present invention provides computational methods for analysis, and synthesis of human motion under partial assist from powered augmentation devices. The algorithms are integrated in a simulation platform to be used as a test-bed for prototyping, simulating, and verifying algorithms to control human motion under artificial control. The analysis and synthesis problems in human motion are formulated as a trajectory tracking control algorithm using inverse and forward models coupled by proportional and derivative feedback terms. A muscle force distribution and capacity module is used to monitor the computed joint torques in order to assess the physiological consequences of the artificial control, and if needed to make modifications. This framework allows us to verify robustness, stability, and performance of the controller, and to be able to quickly change parameters in a simulation environment. We can study many different motions in a simulation environment. Thus, future performance and designs of human augmentation devices can be studied through simulation, without the risk and constraints imposed by hardware implementations of current technology.

[0037] The System Model

[0038] The system (or plant) refers to a dynamic model of the combined musculoskeletal and augmentation device system. The system may be designed having various degrees of complexity, depending on the requirements imposed by the study. Without loss of generality, we consider a simple planar biped system to illustrate the concepts (See FIG. 1). The equations of motion are formulated in such a way to handle three phases of biped motion as shown in FIG. 1. They include single support, double support, and airborn. Let q be the coordinates corresponding to the rotational and translational degrees of freedom.

[0039] where (x₃, y₃) corresponds to the center of mass of the torso and the joint angles are measured clockwise from the vertical.

[0040] The system is actuated by voluntary control from the muscles and artificial control from the augmentation device. The total torque applied at the joints (net joint torque) are the combined torque from the muscles (_m) and the assist actuators (_a)

[0041] Let C(q) represent the foot-floor contact constraints and _LR^T be the vector corresponding to the ground reaction forces under the left and right feet. The equations of motion of the system are given by, $\begin{array}{cc}J\left(q\right)\stackrel{}{q}+B\left(q,\stackrel{.}{q}\right)\stackrel{.}{q}+G\left(q\right)+T_{\mathrm{ad}}=\frac{C^T}{q}{}^{}+D& \left(3\right)\end{array}$

[0042] where J, B, and G correspond to the inertia, coriolis and centrifugal torques, and gravitational terms, respectively. The vector T_ad models the augmentation device dynamics and the constant matrix D characterizes the torque coupling effects at the joints. The matrix D is present because absolute coordinates for the joint angles are used in the formulation of the equation of motion, as opposed to relative coordinates. The ground reaction forces may be expressed as a function of the state and inputs by (Hemami, H., A feedback On-Off Model of Biped Dynamics. IEEE Transactions on Systems, Man, and Cybernetics, Vol. SMC-10, No. 7, July 1980). $\begin{array}{cc}\begin{array}{c}={\left(\frac{C}{q}{J\left(q\right)}^{-1}\frac{C^T}{q}\right)}^{-1}(-\frac{}{q}\left(\frac{C}{q}\stackrel{.}{q}\right)\stackrel{.}{q}+\\ \frac{C}{q}{J\left(q\right)}^{-1}\left(B\left(q,\stackrel{.}{q}\right)\stackrel{.}{q}+G\left(q\right)+T_{\mathrm{ad}}-D\right))\end{array}& \left(4\right)\end{array}$

[0043] FIG. 2 shows a system model description with intermittent contact of left and right feet with the ground. Forward dynamic simulations are performed by computing the induced accelerations {umlaut over (q)} obtained from Equation 3 and Equation 4, using the simulated state variable q and {dot over (q)} which are obtained by numerical integration.

[0044] The Internal (Inverse) Model

[0045] It has been demonstrated that the behavior of the human body when coupled with a novel mechanical system is very similar to the behavior that results when the controller relies on an internal model. One such internal model is thought to be a forward model, a term used to describe the computations involved in predicting sensory consequences of a motor command. There are a number of studies that have suggested that a forward model may be used by the human central nervous system (CNS) to estimate sensory consequences of motor actions (Wolpert, D. M., Miall, R. C., Kerr, G. K., Stein, J. F. Ocular limit cycles induced by delayed retinal feedback. Exp Brain Res., 96: 173-180, 1993; Flanagan. J. R., Wing, A. M. The role of internal models in motion planning and control: evidence from grip force adjustment during movements of hand held loads. J. Neurosci, 17:1519-1528, 1997). This theory is easily understood when considering transmission delays inherent in the sensory-motor loop. Although a forward model is particularly relevant to feedback control of time delayed systems, an inverse model is sometimes considered to predict the motor commands that are appropriate for a desired behavior (Atkeson, C. G. Learning arm kinematics and dynamics. Annu Rev. Neurosci, 12:157-183, 1989; Kawato, M., Adaptation and learning in control of voluntary movement by the central nervous system. Advanced Robotics 3, 229-249, 1989; Shadmehr, R., Leaming virtual equilibrium trajectories for control of a robot arm. Neural Comput, 2:436-477, 1990; Gomi, H., Kawato, M., The cerebellum and vor/okr learning models. Trends Neurosci, 15:445-453, 1992).

[0046] Inverse models are generally not considered for control of time delayed systems since the controller would seem to not have the ability to respond to the error and results in instability. However, it is plausible that local or intrinsic feedback mechanisms in conjunction with an inverse model can function to stabilize a system with latencies. Local feedback with stabilizing characteristics is believed to exist in humans in the form of viscoelastic properties of muscles and spinal reflex loop. The concept of an inverse model is also attractive for analysis problems of biomechanical quantities, whereby internal loads are estimated from kinesiological measurements. The approach adopted here in developing a computational model of human sensory motor control is based on the concept of an inverse model coupled with nonlinear feedback (FIG. 3). This mechanism is compelling from the standpoint of biomechanical analysis of human motion as well as the synthesis of artificial control. Let q_d represent the desired kinematics, obtained from motion capture data. The following control law (D), when applied to the system equations, will result in a simulated response that will track and reproduce the desired kinematics data, $\begin{array}{cc}D{}^{}=J\left(q\right){\stackrel{}{q}}^{*}+B\left(q,\stackrel{.}{q}\right)\stackrel{.}{q}+G\left(q\right)+T_{\mathrm{ad}}-\frac{C^T}{q}{}^{}& \left(5\right)\end{array}$

[0047] where,

[0048] The diagonal matrices K_p and K_v represent the position and velocity feedback gains, respectively. The eigenvalues of the closed loop system are related to the feedback gains by the following,

[0049] A critically damped response (fastest possible non-oscillatory response) to the tracking error can be achieved by specifying the eigenvalues to be equal, real, and negative. The parameter a is constant and set to 0 or 1, depending on the severity of noise in the measurements. If the desired trajectories are obtained from noisy motion capture measurements, it may be appropriate to set a0 and to specify the eigenvalues to be large and negative. This way, tracking is achieved without the need to compute unreliable accelerations from noisy kinematics data.

[0050] Muscle Force and Muscle Capacity

[0051] The muscle force and muscle capacity module should ideally be implemented in the forward path of the closed loop system (as shown in FIG. 5). However, it may also be implemented as a separate module whose output is used for analysis purpose only. In the latter case, the module's inputs would tap into the required variables of the closed loop system, but the module would not alter the closed loop dynamics.

[0052] In either case, a number of different muscle force distribution algorithms may be implemented. The underlying concepts of our choice of muscle force distribution algorithm is presented below.

[0053] The relationship between the net muscular moment _m and the muscle forces F_m is given by, $\begin{array}{cc}D{}_m=-\frac{L^T}{q}{F}_m& \left(9\right)\end{array}$

[0054] where L is the overall length of the muscle actuator, and L^T/q is an (nm) muscle moment arm matrix. Since the number of muscles (m) exceeds the degrees of freedom (n), the computation of the muscle actuator's excitation inputs (and the resulting forces) from an inverse dynamics computation amounts to solving a problem that is inherently ill-posed. Static, nonlinear optimization has been used extensively to predict the individual muscle forces to produce the required torque. There are several compelling reasons for using static optimization to predict the individual muscle forces: first, static, non-linear optimization techniques have well developed theoretical foundations. With the advance of commercial software for solving general, constrained, multi-variable non-linear optimization problems, it is now possible to solve sophisticated problems numerically in relatively short time. Second, the notion that muscle forces are controlled in some way to optimize physiological criteria has great intuitive appeal. It has been shown that for motions like walking, static optimization yields very similar results to dynamic optimization (Anderson, F C and Pandy, M G., Static and Dynamic Optimization Solutions for Gait are practically equivalent, Journal of Biomechanics 34, 2001, 153-161, 2001).

[0055] A muscle force and muscle capacity module takes the computed torques from the inverse model (denoted by D) as inputs and calculates the muscle forces based on a static optimization criterion (module 410 in FIG. 4). While any cost function can be defined in solving the optimization problem, the one used here minimizes the sum of muscle activations squared $\begin{array}{cc}J=\underset{i=1}{\overset{m}{}}{a}_i^2& \left(10\right)\end{array}$

[0056] where m is the number of muscles crossing the joint, a_i is the activation level for muscle i and is constrained to be between 0.01 and 1.0. A muscle force F_i for muscle i can be represented as below; $\begin{array}{cc}{F}_i=a_i{F}_i^0& \left(11\right)\end{array}$

[0057] where $F_i^0$

[0058] represents a maximum force limit for muscle i. A gradient based technique can be used to numerically solve for the muscle activations that minimize the cost function J while satisfying the joint moment equilibrium for all degrees of freedom of interest. The optimization problem can be solved using constrained nonlinear optimization (Sequential Quadratic Programming; AEM Design). Once the muscle activations are obtained, the muscle force can then be determined using the force-length-velocity-activation relationship of muscle (Zajac, F. E. Muscle and tendon: Properties, models, scaling, and application to biomechanics and motor control. Critical Reviews in Biomedical Engineering, 17(4):359-41 1, 1989; Anderson, F C and Pandy, M G., Static and Dynamic Optimization Solutions for Gait are practically equivalent, Journal of Biomechanics 34, 2001, 153-161, 2001 ; Hungspreugs, P., Thelen, D., Dariush, B., Ng-Thow-Hing, V., Muscle Force Distribution Estimation Using Static Optimization Techniques. Technical Report-Honda R&D Americas, April 2001).

[0059] The computed muscle forces are then compared with physiological capacity of the muscle in the muscle capacity module 420. The maximum force limits can be ascertained from the well-studied force-length-velocity relationship of muscle (Zajac 1989, cited above). In addition, the muscle forces with and without the assist torque are compared in order to assess whether the assist torque control has helped (improved efficiency) or hindered the motion. If the assist torque control hinders motion, the muscle forces are adjusted and feasible joint torques are computed (modules 430 and 440 in FIG. 4). A poorly designed assist control would then result in DD, producing a simulated response that would not track the desired response. If the assist torques are well designed, DD and the resulting motion would track the desired motion.

[0060] Augmentation Device Controller

[0061] The inputs to the human augmentation device may include the sensed state variables q_s and/or {dot over (q)}_s, which can be directly measured or estimated. These inputs, denoted by (q_s, {dot over (q)}_s) represent a subset of the total number of state variables (q, {dot over (q)}) in our human model. In addition to the sensed state variables, measurements may also be used as input to the augmentation device controller. The augmentation device controller output represents the assist torque _a, which is then input to the inverse model.

[0062] Different control strategies may be used by the human augmentation device controller. For example, gravity compensation control can be used for tasks requiring an increase in potential energy of the total system (human and exoskeleton). Such tasks would include lifting objects, carrying loads, climbing stairs, rising from a chair, etc. A different control strategy, or hybrid control strategies, may be suitable for other tasks such as walking or running. Here, we will present the gravity compensation control algorithm.

[0063] By using the Lagrangian, we can assess the total potential energy of the musculoskeletal system. Let U denote the total potential energy stored in the system, $\begin{array}{cc}U=\underset{i-1}{\overset{n}{}}{m}_j{g}^T{X}_1& \left(12\right)\end{array}$

[0064] The torque at joint i due to gravity can be computed by taking the partial derivative of U with respect to q_i, $\begin{array}{cc}D{}_{\mathrm{gr}}=\frac{U}{q_i}=\underset{j=1}{\overset{n}{}}{m}_j{g}^T\frac{X_j}{q_i}& \left(13\right)\end{array}$

[0065] where g^T represents the gravitational acceleration vector, and X_j represents the coordinates of the center of mass of segment j. Suppose the knee joint between segment 1 and segment 2 is actuated by an augmentation device and the angle corresponding to q ₂ (represents q_sq) is measurable. The following control law may be used as one algorithm for the augmentation device controller $\begin{array}{cc}D{}_a=D{}_{gr}=\frac{U}{q_2}=\underset{j=1}{\overset{n}{}}{m}_j{g}^T\frac{X_j}{q_2}& \left(4\right)\end{array}$

[0066] Note that the above control algorithm requires the center of mass positions of all the link segments (denoted by X_j). Although X_j can be derived from measurement of joint angles and segment lengths, it may not be feasible to measure all joint angles and all segment lengths. Alternatively, if the vertical component of the ground reaction force under each foot can be measured or estimated, it is possible to derive an iterative ground up gravity compensation algorithm which would eliminate the need for access to center of mass of every segment.

[0067] Integration of Modules

[0068] The block-diagram of the integrated modules as has been presented in the description is shown in FIG. 5. The Augmentation device controller is presumed to have as inputs the sensed states and output the assist torques. The overall framework is very general and enables flexible design of the augmentation device control signals. The details of such designs are easily made by those skilled in the art.

[0069] FIG. 6 shows a flowchart illustrating a simulation process according to an embodiment of the present invention. At step S605, time t is set to 0. At step S610, desired kinematic data for the combined musculoskeletal and augmentation device system are obtained. The desired kinematic data may be obtained from motion capture data.

[0070] At step S612, the simulated kinematic data is fed back to obtain tracking error.

[0071] At step S615, modified accelerations {umlaut over (q)}* are computed using Equation 6.

[0072] At step 617, the sensed kinematic data is fed back.

[0073] At step S620, assist torques D_a are computed using the augmentation device controller 500.

[0074] At step S625, torques D are computed using Equation 5 (inverse model 300).

[0075] At step S630, muscle forces are checked and adjusted to modify the corresponding torques (muscle force and capacity module 400).

[0076] At step S635, the induced accelerations {umlaut over (q)} are computed using Equations 3 and 4 and the simulated kinematic data q and {dot over (q)} are obtained by numerical integration (modules 200, 210 and 220).

[0077] At step S640, time t is incremented and at step S645, whether t is less than t_c or not is determined. If t is less than t_c, the process returns to step S610. If t is equal to or greater than t_c, the process ends.

[0078] It should be noted that the above-mentioned equations, modules or functions can be implemented in any kind of computing devices, including general-purpose computers such as personal computers, work stations and main frame computers, and ASICs (Application Specific Integrated Circuits).

[0079] In an embodiment a general-purpose computer is employed to implement the invention. The general-purpose computer comprises software representing the above-mentioned equations, modules or functions. The software is preferably contained in computer readable mediums. Computer readable mediums include read only memories, random access memories, hard disks, flexible disks, compact disks and so on.

[0080] Simulations

[0081] A very simple simulation of the tracking system is carried out to illustrate some of the concepts proposed in the description.

[0082] A simulation illustrating the tracking characteristics of the proposed method without acceleration estimates of the reference trajectory, is provided. In particular, the double support phase of the biped system during a squatting maneuver was simulated. The results are illustrated in FIGS. 7 to 9.

[0083] In FIG. 7, the desired and simulated joint trajectories illustrate the effectiveness of the tracking procedure. These results were obtained by setting a0, i.e. no acceleration estimates were used as inputs to the inverse model. The corresponding joint torques and ground reaction forces are depicted in FIG. 8 and FIG. 9, respectively.

[0084] It should be noted that those skilled in the art can modify or change the above-mentioned embodiments, without departing from the scope and spirit of the present invention. It should therefore be noted that the disclosed embodiments are not intended to limit the scope of the invention, but only to exemplarily illustrate the invention.