Dynamic interaction-oriented subject's limb time-varying stiffness identification method and device

Abstract

The disclosure provides a dynamic interaction-oriented subject's limb time-varying stiffness identification method and device. The method includes: the combination of subject's limb displacement and measured force data or the combination of angle and measured torque data is collected; based on the time-varying dynamic system constructed based on a second-order impedance model, the linear parameter varying method is utilized to substitute the time-varying impedance parameters and reconstruct the restoring force/torque expression; iterative identification is performed on variable weights, dynamic interaction force/torque, and restoring force/torque by using time-varying dynamic parameters based on the dynamic interaction force/torque expression expanded from basis function; the time-varying stiffness is solved by using variable weights and dynamic interaction force/torque according to expression with substituted the time-varying impedance parameters. The disclosure not only improves the accuracy of the time-varying stiffness identification technology but also expands the application scenarios of the time-varying stiffness identification technology.

Claims

1. A dynamic interaction-oriented subject's limb time-varying stiffness identification method, comprising the following steps: in a random small-perturbation motion environment, a subject keeping a limb posture unchanged, and by changing states of muscles, a combination of a subject's limb displacement and a measured force data or a combination of an angle and a measured torque data being collected by a data acquisition system; based on a restoring force/torque expression with substituted time-varying impedance parameters and an expression of a dynamic interaction force/torque expanded by a basis function, and according to the combination of the subject's limb displacement and the measured force data or the combination of the angle and the measured torque data, iterative identification of the dynamic interaction force/torque and a restoring force/torque being performed by using parameters of time-varying dynamics; a linear coefficient matrix and the dynamic interaction force/torque obtained in the iterative identification of the parameters being utilized to calculate a time-varying stiffness of the subject's limb according to the restoring force/torque expression with the substituted time-varying impedance parameters; wherein the method for obtaining the restoring force/torque expression with the substituted time-varying impedance parameters comprises the following steps: according to time-varying impedance characteristics of the subject's limb, a second-order impedance model being utilized to construct a time-varying dynamic model of the subject's limb; a linear parameter varying method being adopted to construct a time-varying damping parameter of the subject's limb and the time-varying stiffness parameter of the subject's limb in the time-varying dynamic model of the subject's limb as a weighted sum of a function with a scheduling variable as an independent variable, and the restoring force/torque being transformed to obtain the restoring force/torque expression with the substituted time-varying impedance parameters; wherein the dynamic interaction force/torque is the scheduling variable; the method for obtaining the expression of the dynamic interaction force/torque expanded by the basis function is: the dynamic interaction force/torque being expressed by an expansion of the basis function, and a basis function sequence in a perturbation duration interval and weights of the basis function being utilized to represent the dynamic interaction force/torque, wherein the data acquisition system comprises a collaborative robotic arm, reflective markers, an optical motion capture system, a six-dimensional force/torque sensor and a display, wherein a tool end of the collaborative robotic arm is fixedly provided with the reflective markers, and is connected with the six-dimensional force/torque sensor, which is fixedly connected to an endpoint/joint of the limb when in use; the optical motion capture system is configured to collect a three-dimensional position of the reflective markers and obtain the displacement/angle information of the subject's limb; the reflective markers are configured to record the limb posture; the collaborative robotic arm is configured to construct the random small-perturbation motion environment in a three-dimensional space; the six-dimensional force/torque sensor is configured to collect a measured force/torque information between the subject's limb and the tool end of the robotic arm; the display provides the subject with visual feedback of real-time data from the six-dimensional force/torque sensor during dynamic interactions.

2. The dynamic interaction-oriented subject's limb time-varying stiffness identification method according to claim 1, wherein the method performs the iterative identification of the dynamic interaction force/torque and the restoring force/torque by using time-varying dynamic parameters comprises the following steps: (1) collected subject's limb displacement being combined with the measured force data or the angle being combined with the measured torque data, and initialize an estimated value of an initial restoring force/torque as 0; (2) based on the combination of the subject's limb displacement and the measured force data or the combination of the angle and the measured torque data, combined with an estimated value of the restoring force/torque under a current iteration, and according to the fact that a measured force/torque being equal to a sum of the restoring force/torque and the dynamic interaction force/torque, an g of the dynamic interaction force/torque being obtained; (3) according to the expression of the dynamic interaction force/torque expanded by the basis function, an intermediate variable of the dynamic interaction force/torque being utilized to identify the weights of the basis function; (4) an estimated value of the dynamic interaction force/torque under the current iteration being obtained by using the weights of the basis function and the sequence of the basis functions in the perturbation duration interval; (5) according to the fact that the measured force/torque being equal to the sum of the restoring force/torque and the dynamic interaction force/torque, an intermediate variable of the restoring force/torque being obtained; (6) based on the restoring force/torque expression with the substituted time-varying impedance parameters, the intermediate variable of the restoring force/torque being utilized to obtain the linear coefficient matrix; (7) according to a displacement parameter matrix and the linear coefficient matrix obtained under the current iteration, the estimated value of the restoring force/torque under the current iteration being obtained; (8) return to step (2) until an actual number of the iterations reaches a predetermined maximum number of the iterations, and the linear coefficient matrix, the estimated value of the restoring force/torque and the estimated value of the dynamic interaction force/torque obtained in the last iteration being output.

3. The dynamic interaction-oriented subject's limb time-varying stiffness identification method according to claim 1, wherein the method for collecting the combination of the subject's limb displacement and the measured force data or the combination of the angle and the measured torque data is as follows: in the random small-perturbation motion environment, the subject's limb keeping the limb posture unchanged, and an endpoint or a joint of the subject's limbs being connected with a collaborative robotic arm; in the perturbation duration interval, the subject's limb completing at least one dynamic force/torque interaction cycle, and the combination of the subject's limb displacement and the force data or the combination of the angle and the torque data being collected; wherein a dynamic force/torque being a regular interaction force/torque in a manner of, in three directions/rotation axes of X, Y and Z of the collaborative robotic arm, increasing a positive interaction force/torque first and reducing the positive interaction force/torque to 0, and then increasing a negative interaction force/torque and reducing the negative interaction force/torque to 0.

4. The dynamic interaction-oriented subject's limb time-varying stiffness identification method according to claim 1, wherein the restoring force/torque expression with the substituted time-varying impedance parameters is:
f.sub.r(t)=φ(t)βφ(t)=[{umlaut over (x)}(t){dot over (x)}.sub.g(0,t){dot over (x)}.sub.g(1,t) . . . {dot over (x)}.sub.g(P,t)x.sub.g(0,t)x.sub.g(1,t) . . . x.sub.g(P,t)]β=[I b.sub.0b.sub.1. . . b.sub.Pa.sub.0a.sub.1. . . a.sub.P].sup.T wherein I is an inertial parameter of the subject's limb, g.sub.i(f.sub.v(t)) is a function with the dynamic interaction force/torque f.sub.v(t) as the independent variable, a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p are weights of a 0 . . . P-th-order function, P is the maximum order of the function, f.sub.r(t) is the restoring force/torque of the subject's limb, x(t) is a displacement/angle of the subject's limb collected by an optical motion capture system, β is a linear coefficient matrix composing of the inertial parameter I of the subject's limbs and the weights a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p of the 0 . . . P-th-order function, φ(t) is a displacement/angle parameter matrix consisting of variables x.sub.g(i, t) and {dot over (x)}.sub.g(i, t), and a second-order derivative {umlaut over (x)}(t) of the displacement/angle of the subject's limb, x.sub.g(i, t) is an intermediate variable constructed from the function g.sub.i(f.sub.v(t)) and the displacement/angle x(t) of the subject's limb, {dot over (x)}.sub.g(i, t) is an intermediate variable constructed from the function g.sub.i(f.sub.v(t)) and a first-order derivative {dot over (x)}(t) of the displacement/angle of the subject's limb, i is the order of a g.sub.i function, and t is time.

5. The dynamic interaction-oriented subject's limb time-varying stiffness identification method according to claim 1, wherein the expression of the dynamic interaction force/torque expanded by using the basis function is as follows: $f_v\left(t\right)=\underset{j=0}{\overset{Q}{.Math.}}{c}_j{h}_j\left(t\right)$ wherein f.sub.v(t) is the dynamic interaction force, h.sub.j(t) is the basis function sequence in the perturbation duration interval, c.sub.j represents a weight of the j-th-order basis function, Q is the maximum order of the basis function; and t is time.

6. A dynamic interaction-oriented subject's limb time-varying stiffness identification device, comprising: a data acquisition system, a restoring force/torque expression acquisition system, an iterative parameter identification system, and a dynamic interaction force/torque expression acquisition system, wherein the data acquisition system is configured to collect a combination of a subject's limb displacement and a measured force data or a combination of an angle and a measured torque data by changing states of muscles in a random small-perturbation motion environment, and the subject keeps a limb posture unchanged; the restoring force/torque expression acquisition system is configured to reconstruct a restoring force/torque expression by substituting a time-varying impedance parameter of the subject's limb with a linear parameter varying method according to a time-varying dynamic system of the subject's limb; the iterative parameter identification system is configured to, based on the restoring force/torque expression with the substituted time-varying impedance parameters and an expression of a dynamic interaction force/torque expanded by a basis function, and according to subject's limb displacement/angle information and measured force/torque information, perform iterative identification of the dynamic interaction force/torque and the restoring force/torque by using parameters of time-varying dynamics; the dynamic interaction force/torque expression acquisition system is configured to expand the dynamic interaction force/torque using the basis function, and the dynamic interaction force/torque is represented by the basis function sequence in a perturbation duration interval and weights of the basis function, wherein the data acquisition system comprises a collaborative robotic arm, reflective markers, an optical motion capture system, a six-dimensional force/torque sensor and a display, wherein a tool end of the collaborative robotic arm is fixedly provided with the reflective markers, and is connected with the six-dimensional force/torque sensor, which is fixedly connected to an endpoint/joint of the limb when in use; the optical motion capture system is configured to collect a three-dimensional position of the reflective markers and obtain the displacement/angle information of the subject's limb; the reflective markers are configured to record the limb posture; the collaborative robotic arm is configured to construct the random small-perturbation motion environment in a three-dimensional space; the six-dimensional force/torque sensor is configured to collect the measured force/torque information between the subject's limb and the tool end of the robotic arm; the display provides the subject with visual feedback of real-time data from the six-dimensional force/torque sensor during dynamic interactions.

7. The dynamic interaction-oriented subject's limb time-varying stiffness identification device according to claim 6, wherein the method for the restoring force/torque expression acquisition system to obtain the restoring force/torque expression with the substituted time-varying impedance parameters comprises the following steps: according to time-varying impedance characteristics of the subject's limb, a second-order impedance model being utilized to construct a subject's limb time-varying dynamic model; a linear parameter varying method being adopted to construct a time-varying damping parameter of the subject's limb and the time-varying stiffness parameter of the subject's limb in a dynamic model as a weighted sum of a function with a scheduling variable as an independent variable, and the restoring force/torque being transformed to substitute the time-varying impedance parameters; wherein the dynamic interaction force/torque is the scheduling variable.

8. The dynamic interaction-oriented subject's limb time-varying stiffness identification device according to claim 6, wherein the restoring force/torque expression with the substituted time-varying impedance parameters is:
f.sub.r(t)=φ(t)βφ(t)=[{umlaut over (x)}(t){dot over (x)}.sub.g(0,t){dot over (x)}.sub.g(1,t) . . . {dot over (x)}.sub.g(P,t)x.sub.g(0,t)x.sub.g(1,t) . . . x.sub.g(P,t)]β=[I b.sub.0b.sub.1. . . b.sub.Pa.sub.0a.sub.1. . . a.sub.P].sup.T wherein I is an inertial parameter of the subject's limb, g.sub.i(f.sub.v(t)) is a function with the dynamic interaction force/torque f.sub.v(t) as the independent variable, a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p are weights of a 0 . . . P-th-order function, P is the maximum order of the function, f.sub.r(t) is the restoring force/torque of the subject's limb, x(t) is a displacement/angle of the subject's limb collected by an optical motion capture system, β is a linear coefficient matrix composing of the inertial parameter I of the subject's limbs and the weights a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p of the 0 . . . P-th-order function, φ(t) is a displacement parameter matrix consisting of variables x.sub.g(i, t) and {dot over (x)}.sub.g(i, t), and a second-order derivative {umlaut over (x)}(t) of the displacement/angle of the subject's limb, x.sub.g(i, t) is an intermediate variable constructed from the function g.sub.i(f.sub.v(t)) and the displacement/angle x(t) of the subject's limb, {dot over (x)}.sub.g(i, t) is an intermediate variable constructed from the function g.sub.i(f.sub.v(t)) and a first-order derivative {dot over (x)}(t) of the displacement/angle of the subject's limb, i is the order of a g.sub.i function, and t is time; the expression of the dynamic interaction force/torque expanded by using the basis function is as follows: $f_v\left(t\right)=\underset{j=0}{\overset{Q}{.Math.}}{c}_j{h}_j\left(t\right)$ wherein f.sub.v(t) is the dynamic interaction force, h.sub.j(t) is the basis function sequence in the perturbation duration interval, c.sub.j represents a weight of the j-th-order basis function, Q is the maximum order of the basis function; and t is time.

9. The dynamic interaction-oriented subject's limb time-varying stiffness identification device according to claim 6, wherein the method for the iterative parameter identification system to perform iterative identification of the dynamic interaction force/torque and the restoring force/torque comprises the following steps: (1) the collected subject's limb displacement being combined with the measured force data or the angle being combined with the measured torque data, and initialize an estimated value of an initial restoring force/torque as 0; (2) based on the combination of the subject's limb displacement and the measured force data or the combination of the angle and the measured torque data, combined with an estimated value of the restoring force/torque under a current iteration, and according to the fact that a measured force/torque being equal to a sum of the restoring force/torque and the dynamic interaction force/torque, an intermediate variable of the dynamic interaction force/torque being obtained; (3) according to the expression of the dynamic interaction force/torque expanded by the basis function, the intermediate variable of the dynamic interaction force/torque being utilized to identify the weights of the basis function; (4) an estimated value of the dynamic interaction force/torque under the current iteration being obtained by using the weights of the basis function and the sequence of the basis functions in the perturbation duration interval; (5) according to the fact that the measured force/torque being equal to the sum of the restoring force/torque and the dynamic interaction force/torque, an intermediate variable of the restoring force/torque being obtained; (6) based on the restoring force/torque with the substituted time-varying impedance parameters, the intermediate variable of the restoring force/torque being utilized to obtain the linear coefficient matrix; (7) according to a displacement parameter matrix and the linear coefficient matrix obtained under the current iteration, the estimated value of the restoring force/torque under the current iteration being obtained; (8) return to step (2) until an actual number of the iterations reaches a predetermined maximum number of the iterations, and the linear coefficient matrix, the estimated value of the restoring force/torque and the estimated value of the dynamic interaction force/torque obtained in the last iteration being output.

10. The dynamic interaction-oriented subject's limb time-varying stiffness identification method according to claim 2, wherein the method for collecting the combination of the subject's limb displacement and the measured force data or the combination of the angle and the measured torque data is as follows: in the random small-perturbation motion environment, the subject's limb keeping the limb posture unchanged, and an endpoint or a joint of the subject's limbs being connected with a collaborative robotic arm; in the perturbation duration interval, the subject's limb completing at least one dynamic force/torque interaction cycle, and the combination of the subject's limb displacement and the force data or the combination of the angle and the torque data being collected; wherein a dynamic force/torque being a regular interaction force/torque in a manner of, in three directions/rotation axes of X, Y and Z of the collaborative robotic arm, increasing a positive interaction force/torque first and reducing the positive interaction force/torque to 0, and then increasing a negative interaction force/torque and reducing the negative interaction force/torque to 0.

11. The dynamic interaction-oriented subject's limb time-varying stiffness identification device according to claim 7, wherein the restoring force/torque expression with the substituted time-varying impedance parameters is:
f.sub.r(t)=φ(t)βφ(t)=[{umlaut over (x)}(t){dot over (x)}.sub.g(0,t){dot over (x)}.sub.g(1,t) . . . {dot over (x)}.sub.g(P,t)x.sub.g(0,t)x.sub.g(1,t) . . . x.sub.g(P,t)]β=[I b.sub.0b.sub.1. . . b.sub.Pa.sub.0a.sub.1. . . a.sub.P].sup.T wherein I is an inertial parameter of the subject's limb, g.sub.i(f.sub.v(t)) is a function with the dynamic interaction force/torque f.sub.v(t) as the independent variable, a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p are weights of a 0 . . . P-th-order function, P is the maximum order of the function, f.sub.r(t) is the restoring force/torque of the subject's limb, x(t) is a displacement/angle of the subject's limb collected by an optical motion capture system, β is a linear coefficient matrix composing of the inertial parameter I of the subject's limbs and the weights a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p of the 0 . . . P-th-order function, φ(t) is a displacement parameter matrix consisting of variables x.sub.g(i, t) and {dot over (x)}.sub.g(i, t), and a second-order derivative {umlaut over (x)}(t) of the displacement/angle of the subject's limb, x.sub.g(i, t) is an intermediate variable constructed from the function g.sub.i(f.sub.v(t)) and the displacement/angle x(t) of the subject's limb, {dot over (x)}.sub.g(i, t) is an intermediate variable constructed from the function g.sub.i(f.sub.v(t)) and a first-order derivative {dot over (x)}(t) of the displacement/angle of the subject's limb, i is the order of a g.sub.i function, and t is time; the expression of the dynamic interaction force/torque expanded by using the basis function is as follows: $f_v\left(t\right)=\underset{j=0}{\overset{Q}{.Math.}}{c}_j{h}_j\left(t\right)$ wherein f.sub.v(t) is the dynamic interaction force, h.sub.j(t) is the basis function sequence in the perturbation duration interval, c.sub.j represents a weight of the j-th-order basis function, Q is the maximum order of the basis function; and t is time.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a schematic view of a dynamic interaction-oriented subject's limb time-varying stiffness identification method provided by an embodiment of the present disclosure.

(2) FIG. 2 is a schematic view of a perturbation test provided by an embodiment of the present disclosure.

(3) FIG. 3 is a force cycle change diagram in the experiment process provided by an embodiment of the present disclosure.

(4) FIG. 4 is a schematic view of a time-varying stiffness ellipsoid of an upper limb end provided by an embodiment of the present disclosure.

(5) FIG. 5 is a schematic view of a time-varying stiffness curve of an upper limb end provided by an embodiment of the present disclosure.

DESCRIPTION OF THE EMBODIMENTS

(6) In order to make the purpose, technical solutions and advantages of the present disclosure clearer, the present disclosure will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present disclosure, but not to limit the present disclosure.

(7) In one aspect, the present disclosure provides a dynamic interaction-oriented subject's limb time-varying stiffness identification method, including a motion data acquisition method and a time-varying stiffness identification method for subject's limb time-varying stiffness identification.

(8) The process of motion data acquisition method for subject's limb time-varying stiffness identification is as follows: Under the random small-perturbation motion scenario constructed by the collaborative robotic arm, the impedance characteristics of the subject's limb is changed by adjusting the activation the states of muscles, and the physical sensors are used to obtain the subject's limb motion data.

(9) The time-varying stiffness identification methods of subject's limb include linear parameter varying method, basis function expansion method and iterative identification method.

(10) In the method, the linear parameter varying method converts the time-varying impedance parameters of the subject's limb into the form of time-invariant parameters, and then uses the basis function expansion method to describe the dynamic interaction force of the subject's limb. Thereafter, the iterative identification method is utilized to decouple the subject's limb restoring force and interaction force, and finally the accurate identification of the time-varying stiffness of the subject's limbs is achieved.

(11) The time-varying stiffness identification method of the subject's limb provided by the present disclosure can not only obtain the time-varying stiffness of the limb endpoint, such as the upper limb endpoint, but also can calculate the time-varying stiffness of each joint of the subject's limb. In addition, the present disclosure is not limited to using the combination of dynamic interaction force and displacement data of the subject's limb to identify the stiffness of the subject's limb. When the database used for stiffness identification of the subject's limb is the combination of torque and angle data, the method provided by the present disclosure is also applicable, only that force is replaced with torque, and displacement is replaced with angle.

(12) The Details are as Follows:

(13) (1) Motion data collection for time-varying stiffness identification of subject's limb

(14) The disclosure uses the collaborative robotic arm to construct a three-dimensional spatial random perturbation motion environment. The subject's limb and the end of the collaborative robotic arm are firmly connected together. A certain number of reflective markers are fixed on the end of the collaborative robotic arm, and an optical motion capture system is used to collect the three-dimensional positions of the reflective markers to obtain the displacement/angle information of the subject's limb. A certain number of reflective markers will be fixed on the subject's limb and torso in order to record the limb posture.

(15) The end of the robotic arm is connected with a six-dimensional force/torque sensor to collect the force/torque information between the subject's limb and the end of the robotic arm.

(16) A computer monitor is placed directly in front of the subject to display the instantaneous force/torque of the six-dimensional force/torque sensor exerted by the subject limb on the end of the robotic arm during the experiment, providing visual feedback for the subject.

(17) During the acquisition of the time-varying stiffness data of the subject's limbs, under the random small-perturbation motion environment established by the collaborative robotic arm, the limb end or joint of the subject are fixedly connected to the collaborative robotic arm, keeping the limb posture unchanged, and the time-varying stiffness of the subject's limb is adjusted by changing the states of the muscles.

(18) The time-varying property of subject's limb stiffness is induced by the interaction force/torque recorded by the six-dimensional force/torque sensor at the end of the collaborative robotic arm. In each of the three directions (or rotation axes) of X, Y, and Z, the positive interaction force/torque is increased first and then decreased, and then the negative interaction force/torque is increased and then decreased. In this manner, when the interaction force/torque is changed in all the three directions (or rotation axis) a dynamic force/torque interaction cycle is completed.

(19) Within the specified perturbation time, the subject's limb needs to complete at least one dynamic force/torque interaction cycle, and the completion time of each cycle is not limited.

(20) The displacement/angle information and force/torque information of the subject's limb collected according to the above method can be used to identify the time-varying stiffness of the subject's limbs.

(21) (2) Establishment of the time-varying dynamics model of the subject's limbs

(22) In order to characterize the time-varying stiffness of the subject's limbs, a second-order impedance model with a compact expression and capable of describing the time-varying characteristics of the subject's limb impedance is applied to model the subject's limb stiffness. The second-order impedance model is written as:
I{umlaut over (x)}(t)+B(t){dot over (x)}(t)+K(t)x(t)=f.sub.r(t)  (1)

(23) In the expression, I, B(t) and K(t) are collectively called impedance parameters, I is the inertial parameter of the subject's limb; B(t) is the time-varying damping parameter of the subject's limb; K(t) is the time-varying stiffness parameter of the subject's limb; x(t) is the displacement/angle of the subject's limb collected by the optical motion capture system; f.sub.r(t) is the restoring force/torque data of the subject's limb; {dot over (x)}(t) is the first-order derivative of displacement/angle of the subject's limb; {umlaut over (x)}(t) is the second-order derivative of displacement/angle of the subject's limb. It can be seen from the above dynamic model that the displacement/angle x(t) and restoring force/torque f.sub.r(t) of the subject's limb under the small-perturbation environment are the key data for identifying the subject's limb stiffness. Specifically, the displacement/angle x(t) of the subject's limb is generally obtained by an optical motion capture system; the data recorded by the force/torque sensor on the end of the robotic arm during the interaction experiment is the measured force/torque f.sub.m(t) which consists of restoring force/torque f.sub.r(t) and interaction force/torque f.sub.v(t), and the relationship between the measured force/torque f.sub.m(t), the interaction force/torque f.sub.v(t), and the restoring force/torque f.sub.r(t) is expressed by:
f.sub.m(t)=f.sub.r(t)+f.sub.v(t)  (2)

(24) Therefore, decoupling the measured force/torque f.sub.m(t) to obtain the restoring force/torque f.sub.r(t) and dynamic interaction force/torque f.sub.v(t) of the subject's limb during the dynamic interaction is the key to obtain the time-varying stiffness identification of the subject's limb.

(25) (3) Identification of time-varying stiffness of subject's limbs

(26) (3.1) Linear expression of time-varying impedance parameters of subject's limbs:

(27) The idea of the linear parameter varying method is to determine some time-varying variables that directly affect the time-varying characteristics of the identified system through prior knowledge as scheduling variables, and then model the unknown time-varying impedance parameters as weighted sum of a series of functions of the scheduling variable. In dynamic interaction, the dynamic interaction force/torque f.sub.v(t) is the predominant variable reflecting the stiffness characteristics of the subject's limb. Therefore, the dynamic interaction force/torque f.sub.v(t) is set as the scheduling variable, the time-varying stiffness K(t) of the subject's limb and the time-varying damping B(t) of the subject's limb are constructed as weighted sum of a series of functions with the scheduling variable as the independent variable:

(28) $\begin{array}{cc}\begin{array}{cc}K\left(t\right)=\underset{i=0}{\overset{P}{.Math.}}{a}_i{g}_i\left(f_v\left(t\right)\right)& B\left(t\right)=\underset{i=0}{\overset{P}{.Math.}}{b}_i{g}_i\left(f_v\left(t\right)\right)\end{array}& \left(3\right)\end{array}$

(29) In the expression, B(t) is the time-varying damping parameter of the subject's limb, K(t) is the time-varying stiffness parameter of the subject's limb, g.sub.i(f.sub.v(t)) is the function with the dynamic interaction force/torque f.sub.v(t) as the independent variable, a.sub.i and b.sub.i are the weights of the i-th-order function, P is the maximum order of the function g.sub.i(f.sub.v(t)); and t is time.

(30) By substituting formula (3) into formula (1), the following equation can be obtained:

(31) 0$\begin{array}{cc}I\stackrel{.Math.}{x}\left(t\right)+\underset{i=0}{\overset{P}{.Math.}}{b}_i{g}_i\left(f_v\left(t\right)\right)\stackrel{.}{x}\left(t\right)+\underset{i=0}{\overset{P}{.Math.}}{a}_i{g}_i\left(f_v\left(t\right)\right)x\left(t\right)=f_r\left(t\right)& \left(4\right)\end{array}$

(32) In the expression, I is the inertial parameter of the subject's limb, g.sub.i(f.sub.v(t)) is the function with dynamic interaction force/torque f.sub.v(t) as the independent variable, a.sub.i and b.sub.i are the weights of the i-th-order function, P is the maximum order of the function, f.sub.r(t) is the restoring force/torque of the subject's limb, x(t) is the displacement/angle of the subject's limb collected by the optical motion capture system, {dot over (x)}(t) is the first-order derivative of the displacement/angle of the subject's limb, {umlaut over (x)}(t) is the second-order derivative of the displacement/angle of the subject's limb; and t is time.

(33) The variables x.sub.g(i, t) and {dot over (x)}.sub.g(i, t) are constructed as follows
x.sub.g(i,t)=g.sub.i(f.sub.v(t))x(t){dot over (x)}.sub.g(i,t)=g.sub.i(f.sub.v(t)){dot over (x)}(t)  (5)

(34) In the expression, g.sub.i(f.sub.v(t)) is the function with dynamic interaction force/torque f.sub.v (t) as the independent variable, x(t) is the displacement/angle of the subject's limb collected by the optical motion capture system, {dot over (x)}(t) is the first-order derivative of the displacement/angle of the subject's limb, i is the order of the function g.sub.i(f.sub.v(t)), and t is time.

(35) The formula (4) can be transformed into the following form:
f.sub.r(t)=φ(t)βφ(t)=[{umlaut over (x)}(t){dot over (x)}.sub.g(0,t){dot over (x)}.sub.g(1,t) . . . {dot over (x)}.sub.g(P,t)x.sub.g(0,t)x.sub.g(1,t) . . . x.sub.g(P,t)]β=[I b.sub.0b.sub.1. . . b.sub.Pa.sub.0a.sub.1. . . a.sub.P].sup.T  (6)

(36) In the expression, I is the inertial parameter of the subject's limb, g.sub.i(f.sub.v(t)) is the function with dynamic interaction force/torque f.sub.v(t) as the independent variable, a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p are the weights of the 0 . . . P-th-order function, P is the maximum order of the function g.sub.i(f.sub.v(t)), f.sub.r(t) is the restoring force/torque of the subject's limb, x(t) is the displacement of the subject's limb collected by the optical motion capture system, β is a linear coefficient matrix composing of the inertial parameters I of the subject's limbs and the weights a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p of the 0 . . . P-th-order function, φ(t) is the displacement/angle parameter matrix consisting of the variables x.sub.g(i, t), {dot over (x)}.sub.g(i, t), and the second-order derivative {umlaut over (x)}(t) of the displacement/angle of subject's limb, x.sub.g(i, t) is an intermediate variable constructed through g.sub.i(f.sub.v(t)), x(t), {dot over (x)}.sub.g(i, t) is an intermediate variable constructed through g.sub.i(f.sub.v(t)) and {dot over (x)}(t), i is the order of the function g.sub.i(f.sub.v(t)), and t is time.

(37) Through the above process, the identification problem of time-varying parameters is transformed into the identification problem of time-invariant parameters through linear parameter varying method.

(38) (3.2) Basis function expansion of dynamic interaction force/torque f.sub.v(t):

(39) In order to construct the time-invariant expression of the dynamic interaction force/torque f.sub.v(t), the basis function expansion method is utilized to express the dynamic interaction force. The basis function expansion used for modeling the dynamic interaction force/torque f.sub.v(t) is as follows:

(40) $\begin{array}{cc}{f}_v\left(t\right)=\underset{j=0}{\overset{Q}{.Math.}}{c}_j{h}_j\left(t\right)& \left(7\right)\end{array}$

(41) In the expression, f.sub.v(t) is the dynamic interaction force, h.sub.j(t) is the basis function sequence in the perturbation duration interval, c.sub.j represents the weight of the j-th-order basis function, Q is the maximum order of the basis function; and t is time.

(42) Since the basis function h.sub.j(t) is only related to time, when the structure and order of the basis function are determined, the weight value c.sub.0, c.sub.j, . . . , c.sub.Q can be obtained by means of the linear regression method.

(43) (3.3) Iterative identification of parameters for time-varying dynamic systems

(44) Based on the above description of the dynamic interaction force/torque f.sub.v(t) and the restoring force/torque f.sub.r(t) of the subject's limb, these two forces are decoupled from the measured force/torque data f.sub.m(t) through a step-by-step iterative identification method. The steps of the iterative identification method are as follows:

(45) 1: Initialization is performed, f.sub.m(t), x(t), {dot over (x)}(t) and {umlaut over (x)}(t) are known, and let the initial estimated value of f.sub.r(t) as {circumflex over (f)}.sub.r.sup.0(t)=0.

(46) 2: For the n-th iteration, the intermediate variable {circumflex over (f)}.sub.v.sup.temp(t)=f.sub.m(t)−{circumflex over (f)}.sub.r.sup.n=1(t) is calculated, and let

(47) ${\hat{f}}_v^{\mathrm{temp}}\left(t\right)=\stackrel{Q}{\underset{j=0}{.Math.}}{c}_j{h}_j\left(t\right)$
according to formula (7), so as to identify the weight c.sub.j, and then the estimated value {circumflex over (f)}.sub.v.sup.n(t) of f.sub.v(t) is obtained through formula (7).

(48) 3: The intermediate variable {circumflex over (f)}.sub.r.sup.temp(t)=f.sub.m(t)=f.sub.m(t)−{circumflex over (f)}.sub.v.sup.n(t) is calculated, according to formula (5) and formula (6), let {circumflex over (f)}.sub.r.sup.temp(t)=φ(t)β, so as to identify the parameter β, and the estimated value {circumflex over (f)}.sub.r.sup.n(t) of f.sub.r(t) is obtained through formula (6).

(49) 4: Let n=n+1, go back to step 2 until n exceeds the predefined maximum iteration number, then the iteration is terminated.

(50) In the expression, f.sub.m(t) is the measured force/torque, f.sub.v(t) is the dynamic interaction force/torque, f.sub.r(t) is the restoring force/torque, {circumflex over (f)}.sub.r.sup.0(t) is the initial value of the restoring force/torque, {circumflex over (f)}.sub.v.sup.temp(t) is the intermediate variable, {circumflex over (f)}.sub.rn.sup.−1 (t) is the estimated value of the restoring force/torque in the n−1-th iteration, {circumflex over (f)}.sub.r.sup.temp(t) is the intermediate variable, β is the matrix composing of the inertial parameter I of the subject's limb, the weights a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p of the 0 . . . P-th-order function, φ(t) is the matrix composing of the variables x.sub.g(i, t) and {dot over (x)}.sub.g(i, t), and the second-order derivative {umlaut over (x)}(t) of the displacement/angle of the subject's limb, {circumflex over (f)}.sub.m.sup.n(t) is the estimated value of measured force/torque after the n-th iteration, {circumflex over (f)}.sub.v.sup.n(t) is the estimated value of the dynamic interaction force/torque after the n-th iteration, {circumflex over (f)}.sub.r.sup.n(t) is the estimated value of the restoring force/torque after the n-th iteration; i and j are the order of the corresponding functions, and t is time.

(51) After the above iterative identification algorithms are performed, the weight values in formula (3) can be extracted from the parameter β, and the weight value and the value of the dynamic interaction force are substituted into formula (3) to obtain the time-varying stiffness K(t) of the subject's limb during the dynamic interaction process.

(52) The present disclosure also provides a dynamic interaction-oriented subject's limb time-varying stiffness identification device, including: a data acquisition system, a restoring force/torque acquisition system, an iterative parameter identification system, and a dynamic interaction force/torque acquisition system.

(53) (1) The Data Acquisition System

(54) The data acquisition system is configured to collect the combination of the subject's limb displacement and measured force data or the combination of angle and measured torque data under the random small-perturbation motion environment where the subject changes the muscle contraction level but keeps the limb posture fixed, so as to construct a database for time-varying stiffness identification of subject's limbs.

(55) The data acquisition system includes a collaborative robotic arm, reflective markers, an optical motion capture system, a six-dimensional force/torque sensor and a computer monitor.

(56) The end of the collaborative robotic arm is fixedly connected with the six-dimensional force/torque sensor and a handle along the rotary axis, the reflective markers are fixedly pasted on the handle, and the handle is fixedly connected to an endpoint/joint of the limb when in use.

(57) The optical motion capture system is configured to collect the three-dimensional position of the reflective markers and obtain the displacement/angle information of the subject's limb.

(58) Reflective markers are configured to record limb posture.

(59) The collaborative robotic arm is configured to construct a random small-perturbation motion environment in a three-dimensional space.

(60) The six-dimensional force/torque sensor is configured to collect the measured force/torque information between the subject's limb and the end of the robotic arm.

(61) The display provides the subject with visual feedback of real-time data from the six-dimensional force/torque sensor during dynamic interactions.

(62) (2) The Restoring Force/Torque Acquisition System

(63) The restoring force/torque acquisition system is configured to construct the expression of restoring force/torque with substituted time-varying impedance parameter of the subject's limb by using the linear parameter varying method.

(64) The method for the restoring force/torque acquisition system to obtain the expression of the restoring force/torque with substituted time-varying impedance parameters includes the following steps:

(65) According to the time-varying impedance characteristics of the subject's limb, the second-order impedance model is utilized to construct the time-varying dynamic model of the subject's limb.

(66) By using the linear parameter varying method, the time-varying damping parameters of the subject's limbs and the time-varying stiffness parameters of the subject's limbs in the dynamic model are constructed as the weighted sum of functions with the scheduling variable as the independent variable, and the restoring force/torque is transformed to obtain the expression of the restoring force/torque with the substituted time-varying impedance parameters. The dynamic interaction force/torque is the scheduling variable.

(67) Specifically, the method for the restoring force/torque acquisition system to obtain the expression of the restoring force/torque with substituted time-varying parameters includes the following steps:

(68) According to the time-varying impedance characteristics of the subject's limb, the second-order impedance model is used to construct the subject's limb time-varying dynamics, and its expression is:
I{umlaut over (x)}(t)+B(t){dot over (x)}(t)+K(t)x(t)=f.sub.r(t)

(69) In the expression, t is time, I is the inertial parameter of the subject's limb; B(t) is the time-varying damping parameter of the subject's limb; K(t) is the time-varying stiffness parameter of the subject's limb; x(t) is the displacement/angle of the subject's limb collected by the optical motion capture system; f.sub.r(t) is the restoring force/torque data of the subject's limb; {dot over (x)}(t) is the first-order derivative of displacement/angle of the subject's limb; and {umlaut over (x)}(t) is the second-order derivative of displacement/angle of the subject's limb.

(70) The linear parameter varying method is adopted to construct the time-varying damping parameter of the subject's limb and the time-varying stiffness parameter of the subject's limb in the dynamic model as the weighted sum of function with the scheduling variable dynamic interaction force/torque, as the independent variable, and its expression is:

(71) $K\left(t\right)=\underset{i=0}{\overset{P}{.Math.}}{a}_i{g}_i\left(f_v\left(t\right)\right)B\left(t\right)=\underset{i=0}{\overset{P}{.Math.}}{b}_i{g}_i\left(f_v\left(t\right)\right)$

(72) In the expression, B(t) is the time-varying damping parameter of the subject's limb, K(t) is the time-varying stiffness parameter of the subject's limb, g.sub.i(f.sub.v(t)) is the function with the dynamic interaction force/torque f.sub.v(t) as the independent variable, a.sub.i and b.sub.i are the weights of the i-th-order function g.sub.i, P is the maximum order of the function; and t is time.

(73) After substituting the time-varying impedance parameters as above, the time-varying dynamics model of the subject's limb is expressed as:

(74) $I\stackrel{.Math.}{x}\left(t\right)+\underset{i=0}{\overset{P}{.Math.}}{b}_i{g}_i\left(f_v\left(t\right)\right)\dot{x}\left(t\right)+\underset{i=0}{\overset{P}{.Math.}}{a}_i{g}_i\left(f_v\left(t\right)\right)x\left(t\right)=f_r\left(t\right)$

(75) In the expression, t is time; I is the inertial parameter of the subject's limb, g.sub.i(f.sub.v(t)) is the function with dynamic interaction force/torque f.sub.v(t) as the independent variable, a.sub.i and b.sub.i are the weights of the i-th-order function g.sub.i, P is the maximum order of the function, f.sub.r(t) is the restoring force/torque of the subject's limb, x(t) is the displacement/angle of the subject's limb collected by the optical motion capture system, {dot over (x)}(t) is the first-order derivative of the displacement/angle of the subject's limb, and {umlaut over (x)}(t) is the second-order derivative of the displacement/angle of the subject's limb.

(76) The intermediate variables x.sub.g(i, t) and {dot over (x)}.sub.g(i, t) are constructed with expressions as below:
x.sub.g(i,t)=g.sub.i(f.sub.v(t))x(t) {dot over (x)}.sub.g(i,t)=g.sub.i(f.sub.v(t)){dot over (x)}(t)

(77) In the expressions, x.sub.g(i, t) is the intermediate variable constructed by the function g.sub.i(f.sub.v(t)) and the displacement/angle x(t) of the subject's limb, {dot over (x)}.sub.g(i, t) is the intermediate variable constructed by the function g.sub.i(f.sub.v(t)) and the first-order derivative {dot over (x)}(t) of the displacement/angle of the subject's limb, i is the order of the g.sub.i function, and t is time.

(78) Based on the above variables and expressions substituted by time-varying impedance parameters, the restoring force/torque can be written as:
f.sub.r(t)=φ(t)βφ(t)=[{umlaut over (x)}(t){dot over (x)}.sub.g(0,t){dot over (x)}.sub.g(1,t) . . . {dot over (x)}.sub.g(P,t)x.sub.g(0,t)x.sub.g(1,t) . . . x.sub.g(P,t)]β=[I b.sub.0b.sub.1. . . b.sub.Pa.sub.0a.sub.1. . . a.sub.P].sup.T

(79) In the expression, I is the inertial parameter of the subject's limb, g.sub.i(f.sub.v(t)) is the function with dynamic interaction force/torque f.sub.v(t) as the independent variable, a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p are the weights of the 0 . . . P-th-order function, P is the maximum order of the function, f.sub.r(t) is the restoring force/torque of the subject's limb, x(t) is the displacement/angle of the subject's limb collected by the optical motion capture system, β is a linear coefficient matrix composing of the inertial parameter I of the subject's limbs and the weights a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p of the 0 . . . P-th-order function, φ(t) is the displacement/angle parameter matrix consisting of the intermediate variables x.sub.g(i, t), {dot over (x)}.sub.g(i, t), and the second-order derivative {umlaut over (x)}(t) of the displacement/angle of the subject's limb, x.sub.g(i, t) is an intermediate variable constructed from the function g.sub.i(f.sub.v(t)) and the displacement/angle x(t) of the subject's limb, {dot over (x)}.sub.g(i, t) is an intermediate variable constructed from the function g.sub.i(f.sub.v(t)) and the first-order derivative {dot over (x)}(t) of the displacement/angle of the subject's limb, i is the order of the g.sub.i function, and t is time.

(80) (3) The Iterative Parameter Identification System

(81) The iterative parameter identification system is configured to, based on the expression of restoring force/torque with substituted time-varying impedance parameters and the expression of dynamic interaction force/torque expanded by basis functions, and according to the combination of the subject's limb displacement and measured force data or the combination of angle and measured torque data, perform iterative identification to decouple dynamic interaction force/torque and restoring force/torque by using parameters of time-varying dynamics.

(82) The method for the iterative parameter identification system to decouple the dynamic interaction force/torque and the restoring force/torque includes the following steps:

(83) (1) The collected subject's limb displacement is combined with the measured force data or the angle is combined with the measured torque data, and initialize the estimated value of the initial restoring force/torque as 0.

(84) (2) Based on the combination of subject's limb displacement and measured force data or the combination of angle and measured torque data, and the estimated restoring force/torque under the current iteration, according to the fact that the measured force/torque is equal to the sum of restoring force/torque and dynamic interaction force/torque, the intermediate variable of dynamic interaction force/torque is obtained.

(85) (3) According to the expression of the dynamic interaction force/torque expanded by the basis function, the intermediate variable of the dynamic interaction force/torque is utilized to identify the weights of the basis function.

(86) (4) The estimate of the dynamic interaction force/torque under the current iteration is obtained by using the weights of the basis function and the sequence of basis functions in the perturbation duration interval.

(87) (5) According to the fact that the measured force/torque is equal to the sum of the restoring force/torque and the dynamic interaction force/torque, the intermediate variable of the restoring force/torque is obtained.

(88) (6) Based on the expression of restoring force/torque with substituted time-varying impedance parameters, the intermediate variable of restoring force/torque is utilized to obtain a linear coefficient matrix.

(89) (7) According to the displacement parameter matrix and the linear coefficient matrix obtained under the current iteration, the estimate of the restoring force/torque under the current iteration is obtained.

(90) (8) Return to step (2) until the actual number of iterations reaches the predetermined maximum number of iterations, and the linear coefficient matrix, the estimated restoring force/torque and the estimated dynamic interaction force/torque obtained in the last iteration are output.

(91) Specifically, the iterative parameter identification system consists of the following steps:

(92) Step: (1): Initialization is performed, f.sub.m(t), x(t), {dot over (x)}(t) and {umlaut over (x)}(t) are given, and set the initial estimate of f.sub.r(t) as {circumflex over (f)}.sub.r.sup.0(t)=0.

(93) Steps (2), (3), and (4): For the n-th iteration, the intermediate variable {circumflex over (f)}.sub.v.sup.temp(t)=f.sub.m(t)−{circumflex over (f)}.sub.r.sup.n−1(t) is calculated, and set

(94) ${\hat{f}}_v^{\mathrm{temp}}\left(t\right)=\underset{j=0}{\overset{Q}{.Math.}}{c}_j{h}_j\left(t\right),$
so as to identify the weight c.sub.j, and the estimated value {circumflex over (f)}.sub.v.sup.n(t) of obtained f.sub.v(t) is reconstructed.

(95) Steps (5), (6), and (7): The intermediate variable {circumflex over (f)}.sub.r.sup.temp(t)=f.sub.m(t)−{circumflex over (f)}.sub.v.sup.n(t) is calculated, set {circumflex over (f)}.sub.v.sup.temp(t)=φ(t)β, so as to identify the parameter β, and the estimated value {circumflex over (f)}.sub.r.sup.n(t) of obtained f.sub.r(t) is reconstructed.

(96) Step (8): Set n=n+1, go back to step 2 until n exceeds the set maximum number of iterations, then the iteration is terminated.

(97) In the expression, f.sub.m(t) is the measured force/torque, f.sub.v(t) is the dynamic interaction force/torque, f.sub.r(t) is the restoring force/torque, {circumflex over (f)}.sup.0.sub.r(t) is the initial value of the restoring force/torque, {circumflex over (f)}.sub.v.sup.temp(t) is the intermediate variable, {circumflex over (f)}.sub.r.sup.n−1 (t) is the estimated value of the restoring force/torque in the n−1-th iteration, {circumflex over (f)}.sub.r.sup.temp(t) is the intermediate variable, β is the linear coefficient matrix composing of the inertial parameter I.sub.end of the subject's limb, the weights a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p of the 0 . . . P-th-order function, φ(t) is the displacement parameter matrix composing of the variables x.sub.g(i, t) and {dot over (x)}.sub.g(i, t), and the second-order derivative {umlaut over (x)}(t) of the displacement/angle of the subject's limb, {circumflex over (f)}.sub.n.sup.m(t) is the estimated value of measured force/torque of the n-th iteration, {circumflex over (f)}.sub.v.sup.n(t) is the estimated value of the dynamic interaction force/torque of the n-th iteration, {circumflex over (f)}.sub.r.sup.n(t) is the estimated value of the restoring force/torque of the n-th iteration; i and j are the order of the function, and t is time.

(98) (4) The dynamic interaction force/torque acquisition system

(99) The dynamic interaction force/torque acquisition system is configured to expand the dynamic interaction force/torque using the basis function, and the dynamic interaction force/torque is represented by the basis function sequence in the perturbation duration interval and the weights of the basis function.

(100) The step for the dynamic interaction force/torque acquisition system to expand the basis function of dynamic interaction force/torque is as follows: the dynamic interaction force/torque is represented by the basis function sequence and the weights of the basis functions in the perturbation duration interval.

(101) Specifically, the expression of dynamic interaction force/torque expanded by basis functions is:

(102) $f_v\left(t\right)=\underset{j=0}{\overset{Q}{.Math.}}{c}_j{h}_j\left(t\right)$

(103) In the expression, f.sub.v(t) is the dynamic interaction force, h.sub.j(t) is the basis function sequence in the perturbation duration interval, c.sub.j represents the weight of the j-th-order basis function, Q is the maximum order of the basis function, and t is time.

EMBODIMENT

(104) In this embodiment, a six degrees of freedom collaborative robotic arm is used to apply a three-dimensional random perturbation to the end of the upper limbs of the human body, thereby constructing a motion scene of random small-perturbation. The stiffness characteristics of the end of upper limbs are changed by changing the activation state of the muscles, the displacement and restoring force data of the end of upper limbs are collected by physical sensors in the perturbation process, and the three-dimensional time-varying stiffness of the end of upper limbs is calculated according to the provided time-varying stiffness identification method. The overall measurement process is shown in FIG. 1, which includes the following steps:

(105) (1) Collection of motion data of the end of upper limbs

(106) The subject sat on a hard chair, and the handle of the six degrees of freedom collaborative robotic arm was located 390 mm in front of and 100 mm below the rotation midline of the subject's right shoulder, as shown in FIG. 2. After the subject held the handle fixed on the end of the robotic arm, the machinery was activated and the robotic arm performed perturbation. The maximum amplitude of perturbation was set to 10 mm, and the frequency range was set to 0-5 Hz. The setting of parameters may ensure that the perturbation process provides accurate data for the subsequent kinetic identification of ends. The subject's upper limbs maintained a fixed posture, and the activation state of the muscles was adjusted to change the stiffness characteristics of the end of limbs. The force sensor fixed at the end of the robotic arm collected the data of the restoring force and dynamic interaction force generated by the upper limb due to the perturbation. The external optical 3D motion capture system collected displacement data of the end of limbs. During the process of motion data collection, a monitor was placed in front of the subject to display the force/torque change of the six-dimensional force/torque sensor on the end of the robotic arm during the experiment, providing visual feedback for the subject. The subject needed to perform a force cycle with a range of ±30 N in the three directions of X, Y and Z respectively. According to the force feedback information provided by the display, the positive interaction force was first increased and then decreased, and then the negative interaction force was increased then reduced, and the force cycle in three directions was completed in the manner described above, which is a dynamic force interaction cycle. During the perturbation duration interval of 20 seconds, the subject needed to complete at least one cycle of dynamic force transformation, and the completion time of each cycle is not limited. The force cycle of the measured force is shown in FIG. 3.

(107) (2) Establishment of the time-varying dynamics model of the end of upper limbs

(108) In order to characterize the time-varying stiffness characteristics of the end of upper limbs of human body, a second-order impedance model with a simple structure and capable of more comprehensively reflecting the time-varying impedance characteristics of the end of upper limbs of human body is applied to the research of the stiffness identification of the end of upper limbs. The dynamic model of the second-order impedance model is typically modeled as:
I{umlaut over (x)}(t)+B(t){dot over (x)}(t)+K(t)x(t)=f.sub.r(t)  (1)

(109) In the expression, t is time, I is the inertial parameter of the end of upper limbs of human body; B(t) is the time-varying damping parameter of the end of upper limbs of human body; K(t) is the time-varying stiffness parameter of the end of upper limbs of human body; x(t) is the displacement of the end of upper limbs of human body collected by the optical motion capture system; f.sub.r(t) is the restoring force/torque data of the end of upper limbs of human body; {dot over (x)}(t) is the first-order derivative of displacement of the end of upper limbs of human body; and {umlaut over (x)}(t) is the second-order derivative of displacement of the end of upper limbs of human body. It can be seen from the above dynamic model that the displacement x(t) of the end of upper limbs of human body and restoring force f.sub.r(t) of the end of upper limbs under the small-perturbation environment are the key data for identifying the stiffness of the end. Specifically, the displacement x(t) of the end of upper limbs of human body is generally obtained by an optical motion capture system; the restoring force f.sub.r(t) of the end of upper limbs is part of the data recorded by the force/torque sensor on the end of the robotic arm. The data recorded by the force/torque sensor on the end of the robotic arm during the dynamic interaction process is the measured force f.sub.m(t), and the relationship between the measured force f.sub.m(t), the dynamic interaction force f.sub.v(t), and the restoring force f.sub.r(t) is expressed as:
f.sub.m(t)=f.sub.r(t)+f.sub.v(t)  (2)

(110) Therefore, decoupling the measured force f.sub.m(t) to obtain the restoring force f.sub.r(t) and dynamic interaction force f.sub.v(t) of the end of upper limbs of human body during the dynamic interaction process is the key to solving the time-varying stiffness identification of the end of upper limbs of human body.

(111) (3) Calculation of time-varying stiffness of end of upper limbs of human body

(112) (3.1) Linear expression of time-varying impedance parameters of end of upper limbs of human body

(113) The idea of the linear parameter varying method is to determine some time variables that directly affect the time-varying characteristics of the system in the identified system through prior knowledge as scheduling variables, and then model the unknown time-varying parameters as weighted sum of a series of functions of scheduling variables. In the dynamic interaction process, it is difficult to directly obtain the active state of all muscles of the subject's upper limbs, but the joint action of subject's upper limbs allows the end of subject's upper limbs and the external environment to exhibit the dynamic interaction force f.sub.v(t), and therefore the dynamic interaction force f.sub.v(t) directly affects the stiffness characteristics of the system. Based on the above analysis, the dynamic interaction force f.sub.v(t) may be set as the scheduling variable. After the scheduling variable is determined, the time-varying stiffness K.sub.end (t) of the end of subject's upper limbs and the time-varying damping B.sub.end (t) of the end of subject's upper limbs are constructed as weighted sum of a series of functions with the scheduling variable as the independent variable:

(114) $\begin{array}{cc}K\left(t\right)=\underset{i=0}{\overset{P}{.Math.}}{a}_i{g}_i\left(f_v\left(t\right)\right)B\left(t\right)=\underset{i=0}{\overset{P}{.Math.}}{b}_i{g}_i\left(f_v\left(t\right)\right)& \left(3\right)\end{array}$

(115) In the expression, B(t) is the time-varying damping parameter of the end of upper limbs of human body, K(t) is the time-varying stiffness parameter of the end of upper limbs of human body, g.sub.i(f.sub.v(t)) is the function with the dynamic interaction force f.sub.v(t) as the independent variable, a.sub.i and b.sub.i are the weights of the i-th-order function, P is the maximum order of the function; and t is time.

(116) By substituting formula (3) into formula (1), the following may be obtained:

(117) $\begin{array}{cc}{I}_{\mathrm{end}}\stackrel{.Math.}{x}\left(t\right)+\underset{i=0}{\overset{P}{.Math.}}{b}_i{g}_i\left(f_v\left(t\right)\right)\dot{x}\left(t\right)+\underset{i=0}{\overset{P}{.Math.}}{a}_i{g}_i\left(f_v\left(t\right)\right)x\left(t\right)=f_r\left(t\right)& \left(4\right)\end{array}$

(118) In the expression, I.sub.end is the inertial parameter of the end of upper limbs of human body, g.sub.i(f.sub.v(t)) is the function with dynamic interaction force f.sub.v(t) as the independent variable, a.sub.i and b.sub.i are the weights of the i-th-order function, P is the maximum order of the function, f.sub.r(t) is the restoring force of the end of upper limbs of human body, x(t) is the displacement of the end of upper limbs of human body collected by the optical motion capture system, {dot over (x)}(t) is the first-order derivative of the displacement of the end of upper limbs of human body, {umlaut over (x)}(t) is the second-order derivative of the displacement of end of upper limbs of human body; and t is time.

(119) The variables x.sub.g(i, t) and {dot over (x)}.sub.g(i, t) are constructed, and set
x.sub.g(i,t)=g.sub.i(f.sub.v(t))x(t){dot over (x)}.sub.g(i,t)=g.sub.i(f.sub.v(t)){dot over (x)}(t)  (5)

(120) The formula (4) may be transformed into the following form:
f.sub.r(t)=φ(t)βφ(t)=[{umlaut over (x)}(t){dot over (x)}.sub.g(0,t){dot over (x)}.sub.g(1,t) . . . {dot over (x)}.sub.g(P,t)x.sub.g(0,t)x.sub.g(1,t) . . . x.sub.g(P,t)]β=[I b.sub.0b.sub.1. . . b.sub.Pa.sub.0a.sub.1. . . a.sub.P].sup.T

(121) In the expression, I.sub.end is the inertial parameter of the human body, g.sub.i(f.sub.v(t)) is the function with dynamic interaction force f.sub.v(t) as the independent variable, a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p are the weights of the 0 . . . P-th-order function, P is the maximum order of the function, f.sub.r(t) is the restoring force of the human body, x(t) is the displacement of the human body collected by the optical motion capture system, β is a linear coefficient matrix composing of the inertial parameters I of the human body and the weights a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p of the 0 . . . P-th-order function, φ(t) is the displacement parameter matrix consisting of the variables x.sub.g(i, t) and {dot over (x)}g(i, t), and the second-order derivative {umlaut over (x)}(t) of the displacement of human body, x.sub.g(i, t) is an intermediate variable constructed through g.sub.i(f.sub.v(t)) and x(t), {dot over (x)}.sub.g(i, t) is an intermediate variable constructed through g.sub.i(f.sub.v(t)) and {dot over (x)}(t), i is the order of the g.sub.i function, and t is time.

(122) Through the above process, the linear parameter varying method transforms the problem of time-varying system parameters identification into the problem of time-invariant unknown parameters identification. This embodiment chooses to use the 0-3rd order Chebyshev polynomial as the function g.sub.i(f.sub.v(t)) describing the relationship between the time-varying impedance parameters of the end of upper limb and the dynamic interaction force, but the form of the function g.sub.i(f.sub.v(t)) is not limited to the Chebyshev polynomial.

(123) (3.2) Basis function expansion of dynamic interaction force f.sub.v(t):

(124) In order to construct the time-invariant expression of the dynamic interaction force f.sub.v(t), the basis function expansion method is utilized to express the dynamic interaction force.

(125) The basis function expansion for the dynamic interaction force f.sub.v(t) is as follows:

(126) $\begin{array}{cc}{f}_v\left(t\right)=\underset{j=0}{\overset{Q}{.Math.}}{c}_j{h}_j\left(t\right)& \left(7\right)\end{array}$

(127) In the expression, f.sub.v(t) is the dynamic interaction force, h.sub.j(t) is the basis function sequence in the perturbation duration interval, c.sub.j represents the weight of the j-th-order basis function, Q is the maximum order of the basis function; and t is time. Since the dynamic interaction force is characterized in smoothness, the B-spline basis function is selected as the basis function h.sub.j(t) because of being able to well describe the amount of smooth change. Based on the requirement for the identification accuracy of dynamic interaction force, it is necessary to use a basis function with higher order. However, if the order is too high, it is likely to cause an increase in the number of unknown parameters and reduce the overall identification result. In this embodiment, the maximum order of the B-spline basis function is set to 30, and this order setting may substantially ensure that the overall variation regularity of the described variable and the subtle variation in some time intervals may be well estimated.

(128) (3.3) Iterative identification for parameters of time-varying dynamic system for end of upper limbs

(129) Based on the above description of the dynamic interaction force f.sub.v(t) and the restoring force f.sub.r(t) of human body, the two parts of the force are decoupled from the measured force data f.sub.m(t) through a step-by-step iterative identification method. The steps of the iterative identification method are as follows:

(130) 1: Initialization is performed, f.sub.m(t), x(t), {dot over (x)}(t) and {umlaut over (x)}(t) are known, and set the initial estimated value of f.sub.r(t) as {circumflex over (f)}.sub.r.sup.0(t)=0.

(131) 2: For the n-th iteration, the intermediate variable {circumflex over (f)}.sub.v.sup.temp(t)=f.sub.m(t)−{circumflex over (f)}.sub.r.sup.n−1(t) is calculated, and set

(132) 0${\hat{f}}_v^{\mathrm{temp}}\left(t\right)=\underset{j=0}{\overset{Q}{.Math.}}{c}_j{h}_j\left(t\right)$
according to formula (7), so as to identify the weight c.sub.j, and the estimated value {circumflex over (f)}.sub.v.sup.n(t) of obtained f.sub.v(t) is reconstructed.

(133) 3: The intermediate variable {circumflex over (f)}.sub.r.sup.temp(t)=f.sub.m(t)−{circumflex over (f)}.sub.v.sup.n(t) is calculated, according to formula (5) and formula (6), set {circumflex over (f)}.sub.r.sup.temp(t)=φ(t)β, so as to identify the parameter β, and the estimated value {circumflex over (f)}.sub.r.sup.n(t) of obtained f.sub.r(t) is reconstructed.

(134) 4: Set n=n+1, go back to step 2 until n exceeds the set maximum number of iterations, then the iteration is terminated.

(135) In the expression, f.sub.m (t) is the measured force, f.sub.v(t) is the dynamic interaction force, f.sub.r(t) is the restoring force, {circumflex over (f)}.sup.0.sub.r(t) is the initial value of the restoring force, {circumflex over (f)}.sub.v.sup.temp(t) is the intermediate variable, {circumflex over (f)}.sub.r.sup.n−1(t) is the estimated value of the restoring force in the n−1-th iteration, {circumflex over (f)}.sub.r.sup.temp(t) is the intermediate variable, β is the matrix composing of the inertial parameter I.sub.end of the end of upper limbs of human body and the weights a.sub.0 . . . a.sub.p and b.sub.0 . . . b.sub.p of the 0 . . . P-th-order function, φ(t) is the matrix composing of the variables x.sub.g(i, t) and {dot over (x)}.sub.g(i, t), and the second-order derivative {umlaut over (x)}(t) of the displacement of the end of upper limbs of human body, {circumflex over (f)}.sub.m.sup.n(t) is the estimated value of measured force of the n-th iteration, {circumflex over (f)}.sub.v.sup.n(t) is the estimated value of the dynamic interaction force of the n-th iteration, {circumflex over (f)}.sub.r.sup.n(t) is the estimated value of the restoring force of the n-th iteration; i and j are the order of the function, and t is time.

(136) After the above iterative identification algorithms are performed, the parameters in the parameter β, various restoring force f.sub.r(t) and dynamic interaction force f.sub.v(t) are obtained, and the above data is transmitted back to formula (3), so that the time-varying stiffness K.sub.end (t) of the end of upper limbs in the dynamic interaction process may be obtained. By performing singular value decomposition on the obtained stiffness matrix, the schematic view of the stiffness ellipsoid of the end of upper limbs shown in FIG. 4 and the curve diagram K (t) of the stiffness of the end of upper limbs shown in FIG. 5 are obtained.

(137) Compared with the conventional technology, the present disclosure has the following advantages:

(138) The dynamic interaction-oriented subject's limb time-varying stiffness identification method provided by the present disclosure is constructed based on the linear parameter varying method, the basis function expansion method and the iterative parameter identification method of the subject's limb time-varying stiffness identification method. The time-varying stiffness research problem of subject's limbs is transformed into a time-invariant research problem through the linear parameter varying method, thereby providing a new idea for the identification of time-varying stiffness. The provided iterative parameter identification can completely and accurately identify the time-varying stiffness of the subject's limbs by means of data in single perturbation, which significantly reduces complexity of the time-varying stiffness research. More importantly, the disclosure breaks through the limitation that the existing time-varying stiffness identification technology relying on single experimental data cannot be applied to the dynamic interaction with the environment. Generally speaking, the method for identifying the time-varying stiffness of a subject's limb provided by the present disclosure improves the accuracy of the time-varying stiffness identifying technique.

(139) The motion data acquisition method for the time-varying stiffness identification of the subject's limb provided by the present disclosure can uniformly guide the experimental data acquisition of the time-varying stiffness research of the endpoint or joint of the subject's limb. The experimental paradigm is simple, the operability is high, and the stiffness characteristics of the endpoint or joint of the subject's limb can be completely changed. Moreover, the recorded data enables accurate resolution of the time-varying stiffness of the subject's limbs.

(140) The present disclosure utilizes a multi-degree-of-freedom cooperative robotic arm to construct a small-perturbation environment to measure the time-varying stiffness of the subject's limb, thereby saving the complicated process of developing a special-purpose measuring device, reducing the difficulty of measuring the time-varying stiffness of the subject's limb, and expanding the application scenario of using the cooperative robotic arm for stiffness measurement, so that the time-varying stiffness measurement of the upper and lower limbs of the subject can be studied with the assistance of the cooperative robotic arm.

(141) Those skilled in the art can easily understand that the above descriptions are only preferred embodiments of the present disclosure, and are not intended to limit the present disclosure. Any modifications, equivalent replacements and improvements made within the spirit and principles of the present disclosure, etc. should all be included within the protection scope of the present disclosure.