# Identification of Hysteresis in Human Meridian Systems Based on NARMAX Model

Authors: Yonghong Tan; Ruili Dong; Hui Chen; Hong He
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

It has been found that the response of acupuncture point on the human meridian line exhibits nonlinear dynamic behavior when excitation of electroacupuncture is implemented on another meridian point. This nonlinear phenomenon is in fact a hysteretic phenomenon. In order to explore the characteristic of human meridian and finally find a way to improve the treatment of diseases via electro-acupuncture method, it is necessary to identify the model to describe the corresponding dynamic hysteretic phenomenon of human meridian systems stimulated by electric-acupuncture. In this paper, an identification method using nonlinear autoregressive and moving average model with exogenous input (NARMAX) is proposed to model the dynamic hysteresis in human meridian. As the hysteresis is a nonlinear system with multivalued mapping, the traditional NARMAX model is unavailable to it directly. Thus, an expanded input space is constructed to transform the multi-valued mapping of the hysteresis to a one-to-one mapping. Then, the identification method using NARMAX model on the constructed expanded input space is developed. Finally, the proposed method is applied to hysteresis modeling for human meridian systems.


---

## Body

## 1. Introduction

It is known that electro-acupuncture is a modern extension of traditional acupuncture which exists more than three thousand years in Chinese history. Based on the theory of the traditional Chinese medicine, it has been found that the acupuncture points can be grouped in lines called meridian system of the human body. Reference [1] reported that the meridian system has significant effect on human health. References [2, 3] illustrated that the meridian system had architecture with many channels letting electrical signals pass through easily. On the other hand, reference [4] has given a review on electrical properties of acupuncture points and meridians.Nowadays, most of the research results on electro-acupuncture are focusing on the analysis of impedance on acupuncture points [2, 3, 5–8]. Note that the meridian system is, in fact, a network with several lines or channels. In each channel, there are several acupuncture points located along a line. It was found that there were some relations among those points in each channel. Therefore, the analysis just depending on the impedance of a single acupuncture-point would not reflect the main characteristic of the signal transmission in human meridian system. One of the alternatives is to stimulate an acupuncture-point in a channel with electrical signal and measure the corresponding responses of the other acupuncture points in the same channel simultaneously. Thus the corresponding signal transmission performance of the measured channel can be evaluated. In order to properly evaluate the properties of the meridian channels, construction of an accurate model to describe the behavior of the meridian channel is necessary. Moreover, the obtained models can also be used for understanding and improving the acupuncture therapy, acupuncture doctor training, therapy scheme evaluation, and so forth.Some research results have illustrated that the human meridian system is a dynamic system [2, 3, 9]. In this case, the identification of the model to describe the dynamic behavior of the meridian is an efficient way for performance evaluation, simulation, and design of optimal stimulation for disease treatment. Reference [9] developed an auto-regressive and moving average model to describe the human meridian system. It fits the response well when the exciting signal with slow frequency and the input amplitude is rather small. However, when the frequency of the exciting input is higher or the amplitude of the exciting signal is larger, it will illustrate some nonlinear behavior especially hysteretic behavior can be found. Thus, a nonlinear dynamic model should be considered to describe this system.In this paper, the dynamic hysteretic behavior of human meridian system is described. As hysteresis is a nonsmooth nonlinearity with multi-valued mapping, the traditional identification which is used for smooth systems with one-to-one mapping is unavailable directly. In order to transform the multi-valued mapping of hysteresis to a one-to-one mapping, an expanded input space is formed. Then, a nonlinear auto-regressive and moving average model with exogenous input is applied to the identification of the dynamic hysteresis of human meridian. This paper is organized as follows: in Section2, the hysteretic phenomenon in human meridian is discussed. Then, the expanded input space constructed by input and a modified hysteretic operator is presented in Section 3. Based on the expanded input space, the NARMAX model to describe the hysteresis behavior of human meridian is developed in Section 4. The recursive general identification algorithm (RGIA) is used for parameter estimation of the hysteresis NARMAX model. In Section 5, the convergence of the identification algorithm is discussed. After that, in Section 6, the experimental results of the proposed modeling method is illustrated. Finally, the brief summary is presented in Section 7.


## 2. Dynamic Hysteresis in Human Meridian

In this section, we will show the corresponding response of the human meridian when the sinusoidal input with sinusoidal signal with different frequencies is applied to the stimulation of meridian system. In this experiment, the Tai-Yin-Lung channel of meridian is chosen for test. The corresponding experimental configuration of meridian signal measurement is shown in Figure1. The stimulation point is Chize (LU5) while point Shaoshang (LU11) is connected to the ground. The measured meridian point is Lieque (LU7). Figure 2 shows that the hysteretic behavior happened in the response of the measured point Lieque. Note that the corresponding hysteresis behavior varies with different frequencies. It is that the widths of the closed curves are increased as the frequency of input signal increases. Thus, it can be seen that the nonlinearity obviously exists in the response due to the phenomenon of asymmetry and rate-dependence. Moreover, from the experimental results, one can know that the hysteresis in human meridian has the following characteristic, that is,(1)
performance dependent on the frequency change of input;(2)
multi-valued mapping of the output against its input; and(3)
nonsmoothness at the turning points when the input with lower frequency is implemented.Figure 1
Experimental configuration of meridian signal measurement.Figure 2
Hysteresis curves at Lieque point.In this situation, the traditional identification methods are rather difficult to be directly used to model the hysteresis existing in human meridian. The existence of hysteresis in human meridian might be explained as tension and spasm caused by external stimulation.Up till now, there have been some methods proposed for hysteresis modeling, for example, Preisach model [10, 11] and Prandtl-Ishlinskii (PI) model [12]. Those methods used simple backlash operators as the basis functions for modeling. Therefore, lots of operators should be employed in order to obtain accurate models. Although there have been some modified Preisach model [10, 13] as well as modified PI model [14] proposed to describe the hysteresis systems, the structures of those modified schemes are still very complex. In order to simplify the architecture of the model to describe the behavior of hysteresis, references [15] as well as [16] developed the so-called expanded input space based hysteretic models. In the expanded input space, a hysteretic operator which extracted the main movement feature of the hysteresis was introduced as one of the coordinates in the expanded input space. Thus, the multi-valued mapping of hysteresis could be transformed to a one-to-one mapping between the inputs and output. One of the advantages of the expanded input space method is that it can avoid the gradient computation especially at the nonsmooth points of the hysteresis.


## 3. Expanded Input Space

In this section, a hysteretic operator is developed to extract the main feature of hysteresis and to form a coordinate besides the input coordinate in the expanded input space.It is well known that hysteresis can be described by the Preisach model. Hence, the hysteresis is continuous and forms a closed loop in the input-output diagram when the input cycles between two extrema.Then, based on [16], the hysteresis operator f(x) is defined as:
(3.1)f(x)=(1-e|x-xp|)(x-xp)+f(xp),
where x is the current input, f(x) is the current output, xp is the dominant extremum adjacent [16] to the current input x. f(xp) is the output of the operator when the input is xp.Therefore, based on the properties of hysteresis operator (referred to [16]), if it is in discrete-time case, Lemmas 3.1 and 3.2 are composed, respectively, that is,Lemma 3.1.
Letx(t)∈C(R), and C(R) are the sets of continuous functions on R. If there exist two adjacent time instants t1 and t2 (t1≠t2), such that x(t1)=x(t2), but x(t1) and x(t2) are not the local extrema, then f(x(t1))≠f(x(t2)).Lemma 3.2.
If there exist two adjacent time instantst1 and t2 (t1≠t2), such that f(x(t1))-f(x(t2))→0, then x(t1)-x(t2)→0.Remark 3.3.
In Lemma3.1, if x(t1)=x(t2), x(t1) and x(t2) are not the local extrema, f(x(t1))=f(x(t2)), then in this case, the extreme between t1 and t2 does not get sampled in the sampling period. In this situation, the hysteretic operator cannot describe the change tendency of hysteresis. For example, suppose t1 and t2 (t1<t2) are two adjacent time instants, in this case, if x(t1)=x(t2) and f(x(t1))=f(x(t2)), but the output of hysteresis corresponding to one of the input values is in the increase zone while the output of hysteresis corresponding to another input value is in the decrease zone. In this case, the local extreme between those two output values of the hysteresis is obviously missed in the sampling period. Hence, the hysteretic operator cannot describe the change tendency of hysteresis if this situation happens.To handle the extremum-missing problem, a modified scheme of the hysteresis operator is proposed in the following theorem.Theorem 3.4.
Ifx(t1)=x(t2), and x(t1) and x(t2) are not the local dominant extremes and f(x(t1))=f(x(t2)), where t1 and t2 are the adjacent time instants and t1<t2, then, the local dominant extrema located in the segment between points (x(t1),f(x(t1))) and (x(t2),f(x(t2))) cannot be sampled within the sampling period [t1,t2]. In this case, the missed local dominant extreme (xm(tm),f(xm(tm))) can be estimated by:
(3.2)xm(tm)=x(t1)+ε(tm-t1),(1>ε>0)
and the corresponding modified hysteretic operator at estimated local dominant extreme is
(3.3)fm(xm(tm))={(1-e-xm(tm)+xp)(xm(tm)-xp)+fm(xp),x(tm)>x(t1),(1-exm(tm)-xp)(xm(tm)-xp)+fm(xp),x(tm)<x(t1),
where Δt=t2-t1, tm=t1+εΔt for maximum estimate and tm=t1-εΔt for minimum estimate.Proof.
Asx(t1)=x(t2), x(t1) and x(t2) are not the local extrema, whilst x(t1) and x(t2) are located in the increase and decrease zones, respectively, based on (3.2) and tm=t1+εΔt, it leads to:
(3.4)x(t1+εΔt)=x(t1)+ε(t1+εΔt-t1).
Then
(3.5)xm(tm)-x(t1)=x(t1+εΔt)-x(t1)=ε2Δt>0.
Hence, xm(tm)>x(t1). Moreover, x(t1)=x(t2). So, it leads to:
(3.6)xm(tm)>x(t2).
Thus, xm(tm)=x(t1)+ε(tm-t1), (1>ε>0,tm=t1+εΔt) is the local maximum extreme between x(t1) and x(t2). Then, the corresponding output of the hysteretic operator is
(3.7)f(xm(tm))=(1-e-xm(tm)+xp)(xm(tm)-xp)+f(xp).
Moreover, we have
(3.8)f(x(t1))=(1-e-x(t1)+xp)(x(t1)-xp)+f(xp),(3.9)f(x(t2))=(1-ex(t2)-xm(tm))(x(t2)-xm(tm))+f(xm(t)),
respectively. It is obvious that f(x(t1))≠f(x(t2)), and it satisfies the condition given in Lemma 3.1.Similarly, if(xm(tm),f(xm(tm)))is a local dominant minimum point, f(xm) can be described by
(3.10)xm(tm)=x(t1)+ε(tm-t1),(1>ε>0,tm=t1-εΔt),(3.11)f(xm(tm))=(1-exm(tm)-xp)(xm(tm)-xp)+f(xp).
Hence, by combining (3.7) and (3.11), the estimated local dominant extreme can be rewritten as
(3.12)fm(xm(tm))={(1-e-xm(tm)+xp)(xm(tm)-xp)+fm(xp),x(tm)>x(t1)(1-exm(tm)-xp)(xm(tm)-xp)+fm(xp),x(tm)<x(t1).
Then, we have the following theorem.Theorem 3.5.
If hysteresis can be described by Preisach model, there exists a continuous one-to-one mappingΓ: R2→R, such that H[x(t)]=Γ[x(t),fm(x(t))].Proof.
Firstly, it is proved thatΓ is a one-to-one mapping. Let us consider the cases in the following.
Case  1. Assume that x(t) is not the local extremum. In terms of Lemma 3.1, if there exist two different time instants t1 and t2, then
(3.13)(x(t1),fm(x(t1)))≠(x(t2),fm(x(t2))).
Therefore, the coordinate (x(t),fm(x(t))) is uniquely corresponding to hysteresis H[x(t)].
Case  2. Suppose that x(t) is the local extremum. In this case, for two different time instants t1 and t2, there will be
(3.14)(x(t1),fm(x(t1)))=(x(t2),fm(x(t2))).Case  3. Suppose x(t1)=x(t2), fm(x(t1))≠fm(x(t2)), where t1 and t2 are the adjacent time instants and t1<t2, the local extrema located in the segment between points (x(t1),fm(x(t1))) and (x(t2),fm(x(t2))) cannot be obtained within the time period [t1,t2]. In this situation, based on Theorem 3.4, we can use (xm,fm(xm)) between those two points to approximate the local dominant extrema.
According to the property of the Preisach-type hysteresis,H[x(t1)]=H[x(t2)]. Then the coordinate (x(t),fm(x(t))) will be uniquely corresponding to hysteresis H[x(t)].
Combining the above-mentioned three situations, it is obtained thatΓ is a one-to-one mapping. Next, it will be verified that Γ is a continuous mapping.
According to the property of the Preisach-type hysteresis, if(3.15)x(t1)-x(t2)⟶0,
it yields
(3.16)H[x(t1)]-H[x(t2)]→0.
Then, considering Lemma3.2, if
(3.17)fm(x(t1))-fm(x(t2))→0,
it can be deduced that
(3.18)x(t1)-x(t2)→0.
Then we haveH[x(t1)]-H[x(t2)]→0.
Therefore, it can be concluded that there exists a continuous one-to-one mappingΓ: R2→R such thatH[x(t)]=Γ(x(t),fm(x(t))).Remark 3.6.
The theorem stated above shows that the constructed expanded input space consisting of two coordinates, that is,x and fm(x).Remark 3.7.
Based on Theorem3.5, the modified hysteretic operator is the combination of (3.2)–(3.4). We can use the modified hysteretic operator to construct an expanded input space.It has also been proved that the expanded input space is a compact set [16]. Hence, the mapping between the output and the input of the hysteresis on this expanded input space is a one-to-one mapping.


## 4. Nonlinear ARMAX Model for Hysteresis

Based on the above-mentioned expanded input space, the nonlinear auto-regressive and moving average model with exogenous input (NARMAX) [17] is applied to the identification of dynamic hysteresis in human meridian. In this model, the output of the hysteretic operator in the expanded input space is used as the exogenous input. The corresponding hysteretic model can be described by:
(4.1)y^(t)=g(U(t-1),Fm(t-1),Y^(t-1)),
where U(t)=[u(t-1),…,u(t-nu)]T is the input vector, Fm(t-1)=[fm(t-1),…,fm(t-nf)]T is the exogenous input vector which is, in fact, the output vector of the hysteretic operator, Y^(t-1)=[y^(t-1),…,y^(t-ny)]T is the auto-regressive vector, and y^(t) is the output of the NARMAX model; g(·) is the dynamic mapping between the input space and the output of the hysteretic model; nu, nf, and ny are the lags for sequences {u}, {fm}, and {y^}, respectively. For simplification of the representation, (4.1) can also be written as
(4.2)y^(t)≜g(x1(t),…,xp(t)),
where p is the sum of the maximum lag of the model, {x1,…,xn} is the set which consists of the auto-regressive, moving average, and exogenous input items of the model. Hence, the corresponding NARMAX model based on the expanded input space can be described by:
(4.3)y^(t)=∑i1=1pθi1xi1(t)+∑i1=1p∑i2=i1pθi1i2xi1(t)xi2(t)+⋯+∑i1=1p⋯∑iq=iq-1pθi1⋯iqxi1(t)⋯xiq(t),
where θ=[θ1,…,θn]Tis the coefficient vector comprised by n=(p+q)!/(p!q!)-1 coefficients, q is the order of the model. Moreover, in order to simplify the representation of the model, all the variables of the model in the right hand side of the model can be combined into a vector h(t)=[h1(t),…,hn(t)]T constituted by the product combinations of x1(t) to xp(t) from degree 1 to q, that is:
(4.4)h1(t)=x1(t),⋮hp(t)=xp(t),hp+1(t)=x12(t),hp+2(t)=x1(t)x2(t),⋮hn(t)=xpq(t).Thus, (4.3) can be described as a linear regression expression:
(4.5)y^(t)=hT(t)θ.For the determination of the model structure, that is, the estimate of the order and lags of the model from the input and output data, the Akaike’s information criterion (AIC) [18] is applied. As the AIC method needs to know the probability distribution of the data which is usually unknown, we replace the likelihood function to the quadratic model error function. Thus, we have:
(4.6)AICm(q^)=-2logQ(θ^)+2q^,
where θ^ and q^ are the estimate of the coefficient vector θ and the model order, respectively; Q(θ^) represents the quadratic model error function under the condition of θ^, that is,
(4.7)Q(θ^)=∑k=1Ne2(k)2N,
where e and N are the model error and the data number, respectively. By finding the maximum value of AICm(q^). Then, we should form a proper combination of nu, nf, and ny via trial and error to the minimize Q(θ^), and satisfy the constraint nv+nf+ny=p. Finally, the obtained minimized Q(θ^) is used to derive the maximum value of AICm(q^). Thus, the appropriate q^ and the proper combination of nv, nf, and ny can be specified.Then, the corresponding recursive general identification algorithm (RGIA) [19] is applied to the estimation of the parameters of the proposed model. This algorithm is to minimize the quadratic criterion, that is,
(4.8)θ^(t)=argminθ∑k=1Nγ(t)[y(t)-h^T(t)θ^(t-1)]2,
where the data vector h(k) is replaced by the corresponding h^(t), γ(t)>0 is a weighted factor, and γ(k)=∑^-1(k); N is the number of data. The online RGIA is as follows:
(4.9)e^(t)=y(t)-h^T(t)θ^(t-1),(4.10)θ^(t)=θ^(t-1)+K(t)e^(t),(4.11)S(t)=h^T(t)P(t-1)h^(t)+μ(t)Σ(t),(4.12)K(t)=P(t-1)h^(t)S-1(t),(4.13)P(t)={(1μ(t))[I-K(t)h^T(t)]P(t-1)[I-K(t)h^T(t)]T+K(t)Σ(t)KT(t)},(4.14)∑^(t)=∑^(t)+ρ(t)[e^2(t)-∑^(t-1)],(4.15)y^(t)=hT(t)θ^(t-1),
where ρ(k) is the convergent factor, which not only satisfies the conditions given by [20] but also is located within (0, 1), and
(4.16)μ(t)=ρ(t-1)ρ(t)[1-ρ(t)],(4.17)P(0)=λI,0<λ<∞.


## 5. Convergence of the Identification Algorithm

Definition 5.1 (see [20]).
Leth^T(t)[θ^(t-1)-θ]=-ξ(t)e^(t), θ is the true value vector of the model parameters. Moreover, ξ(t) means that the identification is affected by the estimated errors of the internal variables.Assumption 5.2.
The hysteretic system existing in human meridian is BIBO stable, soy(t) is bounded for any bounded input u(t).In the following, the convergence analysis on the recursive identification algorithm in the case of absence of noise will be discussed.Lemma 5.3.
Assume input signalu(t) that is bounded but can fully excite the system, then it satisfies α≤(1/n)∑t=1nh^T(t)h^(t)≤β, where both α and β are bounded positive real numbers.Proof.
Suppose that input signalu(t) is bounded. Considering Assumption 5.2 leads to the hysteresis system of human meridian is BIBO stable. Thus, the elements of the estimated data vector, that is, h^(t) are bounded. That implies that there exist two positive real numbers α and β, and they lead to
(5.1)α≤1n∑t=1nh^T(t)h^(t)≤β.Theorem 5.4.
For the algorithm described by (4.9)–(4.15), which is applied to the parameter estimation of the hysteretic model represented by (4.1)–(4.5), if it matches the following conditions, that is,(1)
the condition of Lemma5.3 is held;(2)
P(k) is a positive definite matrix; moreover,(3)
(5.2)limk→∞λmax[P(k)]→0,(5.3)limk→∞supλmax[P(k)]λmin[P(k)]<∞,(5.4)1-Δ1(k)≤ξ(k)≤1+Δ1(k),then the estimated parameter vector θ^(k) will converge to the true value of the parameter vector θ, as k→∞, where λmax[P(k)] and λmin[P(k)] are the maximum and minimum eigenvalues of P(k), respectively. In (5.4), Δ1(k)is equal to
(5.5)1-h^T(k)P(k-1)h^(k)h^T(k)P(k-1)h^(k)+Σ(k).
The proof of this theorem is given in Appendix.


## 6. Experimental Results

The proposed identification method is applied to modelling of hysteresis inherent in human meridian systems. In this experiment, the sinusoidal signal with attenuating amplitude shown in (6.1) is used to stimulate the meridian point Kongzui on the Taiyin-lung channel but measured the response of another meridian point Chize in the same channel:
(6.1)u=3e-1.75tsin60πt+0.5,
while the input signal
(6.2)u=3e-2.5tsin60πt+0.5
is implemented to excite the same point and obtain the data for model validation test.The corresponding input-output plot of the obtained hysteretic operator mentioned in Section3 is shown in Figure 3, and the hysteresis operator can be used to extract the tendency of hysteresis. Then, the corresponding hysteresis NARMAX model is used to describe the human meridian channel. Moreover, based on the derived input and output data, the model structure is determined by the AIC criterion represented by (4.6) and (4.7). In the experiment, the order-selection procedure is demonstrated in Table 1. It shows that the best result is achieved when nu=1, nf=1, and ny=0, respectively. By randomly setting the initial values of the coefficient vector within [0,1] and λ=106 of the covariance matrix for the RGIA shown in (4.17), the derived NARMAX model validation result for hysteresis in human meridian is shown in Figure 4 while Figure 5 illustrates the corresponding modelling error of the NARMAX hysteretic model.Table 1
Model order selection.
Model order
Value of loss function
n
y
=
2, nu=2, nf=2
0.00035181
n
y
=
2, nu=1, nf=2
9.1954
e
-
006
n
y
=
1, nu=1, nf=2
9.1853
e
-
006
n
y
=
1, nu=1, nf=1
2.5496
e
-
005
n
y
=
2, nu=1, nf=1
8.9733
e
-
006
n
y
=
1, nu=2, nf=1
8.846
e
-
006
n
y
=
0, nu=1, nf=1
8.6718
e
-
006Figure 3
The input-output plot of the hysteretic operator used in expanded input space.Figure 4
Model validation of the NARMAX model for hysteresis in meridian.Figure 5
The modeling error of the NARMAX hysteretic model.For comparison, the ARMAX model on the expanded input space is also used to identify the hysteresis inherent in human meridian. The model structure is chosen asny=0, nu=2, and nf=2, respectively. The recursive least squares (RLS) algorithm is implemented for parameter estimate. Figures 6 and 7 illustrate the model validation result and the modeling error of the ARMAX model, respectively. From the comparison, it is seen that the NARMAX model has achieved better modeling result than that of the ARMAX method.Figure 6
Model validation of the ARMAX model for hysteresis in meridian.Figure 7
The modeling error of the ARMAX hysteretic model.


## 7. Conclusion

In this paper, a nonlinear auto-regressive and moving average model with exogenous input is applied to modeling the hysteretic behavior of human meridian systems. As hysteresis is a nonsmooth nonlinearity with multi-valued mapping, the traditional identification method which is used for smooth systems with one-to-one mapping is unavailable to this case directly. In order to transform the multi-valued mapping of hysteresis to a one-to-one mapping, an expanded input space is formed. Then, a nonlinear auto-regressive and moving average model with exogenous input is applied to the identification of the dynamic hysteresis of human meridian. After that, the RGIA method is used for parameter estimate of the model. The experimental results show that the proposed method is better than the ARMAX modeling method. It provides us with a way to identify the model of the meridian channel in human body so as to construct the potential digital human meridian system for doctor training, evaluation of disease treatment scheme, and human meridian simulations.



---

*Source: 451927-2012-11-01.xml*