Abstract
We introduce an extension of the mass-spring model applied to human tissue that can truly reflect the mechanical properties of a living person's tissue under large deformations under pressure by a massage head. The model is used in the force-position control system and helps us to control the massage force more accurately. Parameters of the model were determined according to the experimental results. Simulation results demonstrated that the model succeeded in modeling living tissue under large deformation, and the force-position control system based on the soft tissue model under large deformation achieved higher accuracy than the control system based on the soft tissue model with constant stiffness.
1. Introduction
Chinese massage has been found to be an efficient way to ease and treat chronic pain [1]. With modernization, massage robots have been used with success in place of traditional massage techniques. Japan's Sanyo Electric Co. first utilized robotic kneading action in 1996 as an alternative for hand massage [2]. Toyohashi University of Technology proposed a four-finger massage hand with force sensor in 2006 [3–6]. Tottori University subsequently introduced a massage chair based on skin resilience in 2007 [7], and in the same year, Japan's Waseda University and Asahi University developed robots for facial massage [8], and Konkuk University applied the robot to back massage [9]. Although these studies discussed the specific aspects of each massage, the entire massage system that integrates massage skills together with the therapy process needs further study. In addition, safety considerations and individualized massage therapy solutions should be incorporated into robotic massage system, so as to make them more adaptable for clinical use.
For these reasons, we developed a robotic massage system acting on acupuncture points based on traditional Chinese massage therapy theory, with reproduced human manipulation technics based on the experience of experts [10]. The kinematic and force features of key massage technics, such as thumb kneading, pressing, rolling, vibrating, and pinching, were summarized by analyzing the massage processes of massage therapists, and a mathematical model for robotic massage was established. A pain threshold value was introduced to individualize therapy schemes, and a force-position control method based on the pain threshold was presented. VAS (Visual Analogue Scale) tests for lumbar muscle strain were carried out using the massage robot to validate the treatment effect.
The robot end-effector was contacted with subjects' soft tissue when working, and the control system was affected, which allowed the proposal of a system based on soft tissue [11]. However, the system was not stable across individuals because the soft tissue model was simplified as a spring-damping system with constant stiffness, which is only applied to people under small deformations. To achieve more precise control of massage strength, soft tissue should be modeled under large deformations. Therefore, this paper presents a method to model living tissue under large deformations and a corresponding control system.
2. System Overview
The system is composed of a massage bed, a positioning platform, a massage end-effector, a visual tracking system, and auxiliary devices, as shown in Figure 1. The massage bed is the structural base of the whole system, and the patient lies on it during robotic massage. Auxiliary devices, such as audio players, physiological information-collecting devices, and pain threshold acquisition devices, are embedded in or attached to the bed [10].

System configuration of the massage robot.
The positioning platform is mounted on the massage bed and is designed to position the end-effector on the massage spot and fulfill manipulation technics when needed. As the positioning platform runs over the patient's back, it consists of two assemblies, that is, a gantry like a moving frame for translational positioning in three perpendicular directions and an arm with four rotational joints so that the position and orientation of the end-effector could be adjusted for different techniques at specific acupuncture points. The maximum workspace of the platform is 180 × 55 × 60 cm, which allows the end-effector to reach every part of patient's body to perform the massage. However, it is difficult to integrate all the different manipulation techniques for various therapy functions, such as rolling, pressing, vibrating, thumb kneading, and pinching into one mechanical manipulator. Therefore, the end-effector is designed to be composed of two parts, that is, a parallel mechanism with five degrees of freedom (DOFs) to realize rolling, pressing, and vibrating techniques and a humanoid massage hand with one DOF for thumb kneading and pinching techniques. According to the specific requirements for different techniques, the massage techniques are reproduced by the end-effector. The visual-tracking device of a stereo camera is installed on an elongated arm extending from the positioning platform and moves along with the platform. Markers are attached to both the patient body and the robot body, so the tracking device can calculate the target position relative to the robot coordinate frame. With its scope of view covering the workspace of the end-effector, the tracking device performs real-time monitoring in case the patient moves; thus, it can ensure safety and precision of the massage operation. The tracking error of the massage robot's vision servosystem is less than 2.5 mm.
An audio device is embedded into the massage bed to help the patient relax, which is advantageous for augmenting the treatment effect according to traditional Chinese massage therapy theory. Physiological information-collecting devices are also provided as an input to the robotic control system for adaptive control according to the patient's current physiological status. The pain threshold acquisition device is a digital glove equipped with force sensors. Before beginning massage, physicians use this glove to press pain points on the patient to determine their pain threshold value, which is sent to the control system to customize the massage force for the specific patient.
The main control strategy of the system is force-position control; it has two control modes in the massage operation process. The position-tracking control is used in free space before the end-effector contacts with the human body. The other is force-tracking control used in restricted space after the end-effector contacts the human body. The two control modes can be shifted according to the contacting force magnitude. Section 4.2.2 describes the errors of force-position control, which are less than 0.05 N and 0.1 mm.
Moreover, safety elements were designed under the consideration that the robot would be in close contact with the patient. For instance, motion DOF redundancy is provided in the positioning platform to avoid interference of the robot body with patients. There are seven DOFs in the positioning platform, and motion sequences are configured to proceed from one axis to another; the vertical prismatic joint is driven after all the other joints move.
3. Material and Methods
3.1. Modeling Living Tissue under Large Deformation Using an Index Model
3.1.1. Experimental Setup
To establish an accurate human soft tissue model, we needed to perform experiments on the soft tissues of living subjects. Properties of the adipose tissue (mainly fat) have been characterized using compression tests of human breast [12–15], heel [16, 17], and forearm, [18], as well as pig tissue [19, 20]. In this paper, we designed an experimental setup to measure the displacement of living tissue when different forces compress a person's arm. Because our experiment was associated with the Chinese massage robot project, study subject safety was assured. There are very few data available in the literature regarding the properties of living human tissue. Furthermore, precise steps were taken to adjust the parameters of the simulation model for further clinical applications.
The experiment setup for measuring the mechanical properties of living tissue is shown in Figure 2. A 6 mm diameter tip was mounted on a force sensor. The velocity was controlled by a servomotor with a large range of speed. The external force on human tissue was obtained by the force sensor, and the three-dimensional coordinates of the marker on human's arm were later written to a file when the tracker communicated to a PC by A/D. The data of the force and the vertical movement of the tip were used to generate the force-displacement curve.

Experimental setup for measuring mechanical properties of living tissue.
3.1.2. Experiment Results
To obtain the force-displacement curve, the tip was moved at a constant velocity to compress the person's arm until they felt pain. At the same time, the force applied to the person's arm and the tip displacement were measured and recorded at a certain frequency. Five different velocities were considered: 5, 10, 15, 20, and 30 mm/s. Five independent environments were carried out for each velocity. Figure 3 shows the displacement-force relationship curves for one person at five different velocities. The figure illustrates that the displacement-force relationship curves were very similar, so we used the velocity of 30 mm/s to study the human soft tissue model.

Displacement-force relationship curve.
3.1.3. Index Model Method
(1) Mass-Spring Model. Mass-spring system is based on physical law. It has been widely and efficiently used in modeling various deformed objects. In this system, the modeled object is discretized into cubes, each containing four tetrahedrons. Each point is in 4 cubes and has 18 points surrounding it [21].
Usually, the spring's elasticity is linear. In the dynamic system, the motion of each single quality key point follows Newton's second law [22], which is expressed by
where m
i
is the mass of point i and x
i
is a three-dimensional vector presenting the position of point i.
The motion equation of the entire system is obtained by integrating all the motion equations of each point in the grid structure. By summing N position vectors of each independent quality point to a 3N-dimensional vector, we obtain an equation that takes the form of
where M, C, K are all tensors of size 3N × 3N, representing the mass, resistance, and the elastic coefficient. x is a 3N-dimensional position vector. Although these matrices are very large, they are also sparse. M and C are diagonal matrices, and K is a banded matrix because the spring tension depends only on the distance between the adjacent quality points. f is a 3N-dimensional vector, representing the total external forces applied to all the quality points. Equation (2) can be transformed into two first-order differential equations expressed by
where v is the velocity vector of the entire quality point system. Then the x and v function expressions of time t can be solved by using different integral methods.
To save computation time, we chose an explicit scheme based on the Euler method to solve the dynamic system. Vertex position at time t + h is derived from the position at time t using the following equations [21]:
where a i (t), v i (t), and x i (t) are the acceleration, velocity, and the position of point i; F i is the external force of point i; K i, j is the elastic modulus of the spring connecting points i and j; L i, j is the free length of the spring; and γ i and m i are the damping coefficient and the mass of point I, respectively.
The convergence condition is when acceleration and velocity are the lowest. In practice, the system achieves balance when the acceleration and velocity are less than specific given values.
(2) Index Modeling Method. Previous models were based on the assumption of linear elasticity theory. However, the results of our experiment (Section 3.1.2) showed that the deformation of live persons' soft tissues was not linearly elastic. To establish the most suitable model for modeling living tissue under large deformations, the expression of the elastic coefficient in formula (2) must be determined. Because human tissue has viscoelasticity, the deformation should have an exponential delay (
where a and b are the parameters to be determined and x is the deformation displacement according to the experiment results. Then, just as in the mass-spring model, the system can be solved using the Euler method.
In order to find the k expression with regard to deformation displacements, we used the Euler method to calculate k value corresponding to different forces. As shown in Figure 4, the calculating steps for each force are as follows.
An initial force was given.
The parameters of m i , γ i , and h were given, and the k value was estimated.
Take the parameters above into the formula and use the Euler method to simulate.
Compare the simulation and experimental results. If the difference is smaller than a certain value, then output k value or adjust k value and simulation from (ii) again until a suitable k value is obtained.
Change the value of force and jump to (ii).

k value calculation.
3.2. Control System Based on the Soft Tissue under Large Deformation
A control system based on the soft tissue model under large deformation mentioned in Section 3.1 is proposed. The soft tissue is simplified as a spring-damping system, which is expressed by
where F(t) is the contact force between the end-effector of the massage robot and the soft tissue, x(t) is the position of the end-effector, K is the stiffness, and B is the damping coefficient.
Calculating the Laplace transform of (6) yields the following equation:
The values of K and B can be determined depending on the simulation results in Section 3.1, and they can be expressed by
The relationship between the motor output torque and the voltage can be expressed by
In the system, the force controller is controlled by Pf and If, while the displacement controller is controlled by Pd, Id, and Dd.
The Z-axis motor control system was modeled using the Simulink toolbox in Matlab, as expressed in Figure 5.

Z-axis motor control system.
4. Results and Discussion
4.1. Index Modeling Results
Only some types of human tissue are deformed when force is applied to the human arm. So 51 × 51 × 40 vertices and 97,500 cubes were used in the system for the test. The sides of all these cubes are 1 mm. To ensure that the simulation results are convergent, after a few tests, the parameters were given as follows: m i = 0.001 g (according to human density), γ i = 0.0005 · 10−3 N·S/m, h = 0.01 s, and the two convergence thresholds were both 0.0005.
Figure 6 shows the comparison between simulation and experimental results. The black curve in Figure 4 plots the simulated deformation displacement under different forces obtained after parameter fitting, the red curve is the experiment measurement used to calculate the elastic coefficient, and the gray curves are the experiment measurements of four different people at the same speed of 30 mm/s. The simulation shows good correspondence to the red experimental curve, and the trend of the experiment curves similar to all four subjects.

Comparison between simulation and experiment results. Simulation results are in black, the referred experiment results are in red, and the other four experiment measurements at the same speed of 30 mm/s are in gray.
Figure 7 shows the comparison between the fitting and calculating results of the elastic coefficient. The black curve in Figure 5 shows the calculating results of k value corresponding to deformation displacements, while the blue curve shows the fitting curves using index function. The data show that the curve using index function fits closely to the calculating data.

Comparison between fitting and calculating results. Calculating results of k value are in black, and index function fitting results are in blue.
The relationship between force and displacement using the constant stiffness model when k = 0.1 and the index model mentioned above were calculated. The simulation comparison results are shown in Table 1; S1 and S2 are the simulation results using the index model and the constant stiffness model, and the S3 row shows the measurement results. The error comparison results are shown in Table 2, of which Error 1 and Error 2 are the errors between the measurement results and the simulation results using the index model and the constant stiffness model, respectively.
Simulation comparison results.
Error comparison results.
From the above tables, we can see that the maximum error between simulation results and measurement results can be decreased from 4.8856 to 0.4894, and the accuracy can be improved from 57.5% to 4.6%. These findings show that the index model succeeded in more accurately expressing the mechanical property of living tissue under large deformation.
4.2. Control System Results
4.2.1. Force Simulation Results
By imputing the simulation results into (8), we can get B = 1.3 NS/mm, and the K values under different forces are shown in Table 3.
K values under different forces.
Figure 8 shows the force simulation results. From the figure, we can see that the both of the errors are less than 0.1% and the system has good stability.

Force simulation results.
4.2.2. Force Measurement Results
Using the experimental setup mentioned in Section 3.1.1, different force values were input to the massage robot, and we compared the robot's given forces (F1) and actual forces (F2). The results are shown in Table 4.
Force data accuracy.
From Table 4, we can see that the actual force is quite close to the given force. The error was less than 5%, which is more accurate than the previously reported value of 10% [11].
Above all, the force-position control system based on soft tissue under large deformation succeeded in accurately controlling the massage force imposed on people.
5. Conclusion
We designed an experiment setup to measure the displacement of living tissue when different forces compressed a person's arm. We introduced an index model that can accurately reflect the biomechanical characteristics of living tissue under large deformation.
The model allowed us to improve the force-position control system using the simulation results of K and B. The experimental results demonstrate that the control system has succeeded in following the tracks of the given force, and the error is less than the control system based on soft tissue with constant stiffness. Moreover, because the K and B were obtained from the proposed soft tissue model when different forces were imposed, the control system is more stable for making estimates in different people.
Footnotes
Acknowledgment
This project was supported by the National Support Project Chinese Massage Robot.
