CN101816822B - Setting method of functional electrical stimulation PID (Proportion Integration Differentiation) parameter double source characteristic fusion particle swarm - Google Patents

Abstract

The invention relates to the field of instruments for extremity rehabilitation by utilizing electric pulse stimulation and provides a setting method of a functional electrical stimulation PID (Proportion Integration Differentiation) parameter double source characteristic fusion particle swarm. With the setting method, the current strength of an FES (Functional Electrical Stimulation) system can be accurately and stably controlled in real time and the accuracy and the stability of the FES system can be efficiently enhanced. The invention adopts the technical scheme that: firstly, a knee joint angle is predicted by utilizing a handle reaction vector (HRV) in the walk aid process; and secondly, a proportion calculus PID parameter is set by utilizing a chaos particle swarm algorithm, the FES current level strength is regulated in real time, and finally self-adaptive online setting of the proportion calculus PID parameter is realized, and the invention is also used for a functional electrical stimulation FES system. The invention is mainly applied to setting the PID parameter in functional electrical stimulation.

Description

Functional electrostimulation pid parameter double source characteristic fusion particle swarm setting method

Technical field

The present invention relates to carry out the instrument field of limb rehabilitating, especially the double source Feature Fusion chaos particle swarm setting method of pid parameter in the functional electrostimulation with electric pulse stimulation.

Background technology

(Functional Electrical Stimulation is to stimulate limb motion muscle group and peripheral nerve thereof through current pulse sequence FES) to functional electrostimulation, recovers or rebuild the technology of the componental movement function of paralytic patient effectively.According to statistics, because the spinal cord regeneration ability is faint, to the spinal cord injury paralysed patient, the effective treatment method that can directly repair damage is not arranged as yet at present, implementing function rehabilitation training is effective measures.Spinal cord injury paralysed patient number increases year by year, and function rehabilitation training is a technology of demanding demand urgently.The sixties in 20th century, Liberson successfully utilizes the electro photoluminescence nervus peronaeus to correct the gait of hemiplegic patient's drop foot first, has started the new way that functional electrostimulation is used to move and Sensory rehabilitation is treated.At present, FES has become the componental movement function of recovering or rebuilding paralytic patient, is important rehabilitation means.Yet how accurate triggering sequential and the pulse current intensity of controlling FES can accurately be accomplished the key problem in technology that the intended function action is still FES with assurance electro photoluminescence action effect.According to statistics; The mode of the triggering of FES control is at present studied still few; And according to action effect and predetermined action deviation; Come adjustment FES stimulus intensity and time sequence parameter automatically with closed-loop control, thereby improved the accuracy and the stability of FES system greatly, but present effectively control method is still among exploring.

Handle retroaction vector (handle reactions vector; HRV) be that in fact the effectiveness that walker offers the patient can be divided into clear and definite independently 3 parts in the process of standing and walking under helping according to walker: sagittal trying hard to recommend into, about to dynamic balance and the power support of upward and downward; This also can be regarded as the patient in fact and keeps the new ideas that the required to external world additional mechanics demand of self normal stand walking proposes; Promptly be that the patient is reduced to centre-point load to the effect of walker is synthetic in the walking process of standing, represent with two mechanics vectors at handle mid point xsect centre of form place respectively, as shown in Figure 1; Vector is at x; Y, the size of making a concerted effort of the durection component on the z axle can characterize the patient respectively and tries hard to recommend into dynamic balance and power support level by what walker obtained.Wherein, the x axle forward that sets of definition coordinate system is patient's dextrad, and y axle forward be patient's a forward direction, z axle forward be the patient on to.Like this, the defined formula of HRV also can be written as:

[HRV]＝[HRV ₁，HRV _r] ^T＝[F _lx，F _ly，F _lz，F _rx，F _ry，F _rz] ^T (1)

At present, the situation when HRV is widely used in supervision patient walks in the electro photoluminescence process prevents that then patient from falling down, and causes the secondary injury.This patent proposes to utilize this parameter prediction knee joint angle, the accurate then levels of current intensity of controlling the FES system, and assurance electro photoluminescence action effect can accurately be accomplished the intended function action, and prevents muscular fatigue.

Ratio infinitesimal analysis (proportional-integral-differential; PID) be a kind of very practical feedback regulation algorithm; It detects according to system or the operation deviation; Proportion of utilization, integration, the required regulated quantity of acquisition of differentiating are widely used in engineering practice so that system is carried out FEEDBACK CONTROL because of it is easy to operate.Especially indeterminate or when being difficult to timely on-line determination, safe closed-loop control can be adopted the PID setting algorithm when the controlled system characterisitic parameter.In the face of the complicacy and the time variation operating environment of muscle, because good stability, the reliable operation of PID have still obtained in the functional electrostimulation field using widely at present.The PID core technology is accurate confirm wherein ratio, integration, differential coefficient, especially in the FES field, system stability is required very strictness, so select particularly important to pid parameter.PID control will obtain controls effect preferably, must adjust ratio, integration and three kinds of control actions of differential, forms in the controlled quentity controlled variable not only to cooperatively interact but also the relation of mutual restriction.

Summary of the invention

For overcoming the deficiency of prior art; The double source Feature Fusion chaos particle swarm setting method of pid parameter in a kind of functional electrostimulation is provided; Can accurately stablize and control systematically strength of current of FES in real time; Improve FES system accuracy and stability effectively, and obtain considerable social benefit and economic benefit.For achieving the above object, the technical scheme that the present invention adopts is: the double source Feature Fusion chaos particle swarm setting method of pid parameter in the functional electrostimulation comprises:

At first, utilize the handle retroaction vector HRV forecasting knee joint angle of walk help process;

Secondly, utilize the chaos particle swarm optimization ratio infinitesimal analysis pid parameter of adjusting, real-time monitoring FES levels of current intensity, the flow process of adjusting is: at first according to three decision variable K of ratio infinitesimal analysis PID _p, K _iAnd K _dThe bound of span; Confirm parameters such as population population size, search volume dimension, and the speed and the position of initialization particle colony, the fitness value that calculates each particle in the population through the corresponding relation of actual joint angles and muscle model output joint angles as appropriate evaluation function utilized then; And its fitness and optimum position fitness value itself made comparisons; And with it as the particle typical value, then in adjustment particle's velocity and other parameters, change the optimum position of particle; Till stable, calculate the K that final best position promptly gets ratio infinitesimal analysis PID _p, K _iAnd K _dThree coefficients; Computing system output yout under the new ratio infinitesimal analysis PID coefficient and with the deviation of muscle model output joint angles after get into the self study and the weighting coefficient self-adjusting of next step neural network again; This process repeatedly; The self-adaptation on-line tuning of final realization ratio infinitesimal analysis pid control parameter, and be used for functional electrostimulation FES system.

Said muscle model output joint angles is the method that adopts PLS, that is:

Be provided with m HRV variable HRV1 ..., HRVm, p M variable, M1 ..., Mp, common i (i=1 ..., the n) data set of individual observed reading, T, U are respectively the composition that from HRV variable and M variable, extracts, and are called the offset minimum binary factor,

Concentrate the linear combination of extracting first couple of composition T1, U1 to be from original variable:

T ₁＝ω ₁₁HRV ₁+…+ω _1mHRV _m＝ω ₁′HRV (4)

U ₁＝v ₁₁M ₁+…+v _1pM _p＝v ₁′M (5)

ω wherein ₁=(ω ₁₁..., ω _1m) ' be model effect weight, v ₁=(v ₁₁..., v _1p) ' be M variable weight is converted into the requirement of said extracted first composition and asks constrained extremal problem:

Wherein t1, u1 are the score vector of first pair of composition of being tried to achieve by sample, and HRV0, M0 are initializaing variable, utilize method of Lagrange multipliers, and the problems referred to above are converted into asks vector of unit length ω ₁And v ₁, make θ ₁=ω ₁' HRV ₀' M ₀V ₁Maximum is promptly asked matrix H RV ₀' M ₀M ₀' HRV ₀Eigenwert and proper vector, its eigenvalue of maximum is θ ₁ ², corresponding unit character vector is exactly the ω that separates that is asked ₁, and v ₁By formula

Obtain;

Next sets up the equation of initializaing variable to T1

$\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="0"\; label="HRV"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">HRV</figure-callout>}}_0={\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1{\mathrm{\&alpha;}}_1^{\&prime;}+E_1\\ {\mathrm{<figure-callout\; id="0"\; label="M"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">M</figure-callout>}}_0={\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1{\mathrm{\&beta;}}_1^{\&prime;}+F_1\end{array}\right.---\left(7\right)$

T wherein ₁Meaning is the same, α ₁'=(α ₁₁..., α _1m), β ₁'=(β ₁₁..., β _1p) be the parameter vector when only a M measures t1, E1, F1 are respectively n * m and n * p residual error battle array, can try to achieve coefficient vector α according to common least square method ₁And β ₁, α wherein ₁Become model effect load capacity;

Can not reach the precision of regression model like first composition that extracts, utilization residual error battle array E1, F1 replace X0, Y0, repeat to extract composition, and the like, supposing finally to have extracted r composition, HRV0, M0 to the regression equation of r composition are:

$\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="0"\; label="HRV"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">HRV</figure-callout>}}_0={\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1{\mathrm{\&alpha;}}_1^{\&prime;}+...+t_r{\mathrm{\&alpha;}}_r^{\&prime;}+{\mathrm{<figure-callout\; id="0"\; label="E\; r\; M"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">E</figure-callout>}}_r& \\ M_{}{}_0={\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1{\mathrm{\&beta;}}_1^{\&prime;}+...+t_r{{\mathrm{\&beta;}}_r}^{\&prime;}+F_r\end{array}\right.---\left(8\right)$

The first step analyze extract in the gained HRV amount composition Tk (k=1 ..., r) regression equation that the M amount is set up r composition, i.e. t are brought in linear combination into _r=ω _K1HRV ₁+ ... + ω _KmHRV _mSubstitution M _j=t ₁β _1j+ ... + t _rβ _Rj(j=1 ..., p), promptly get the regression equation M of standardized variable _j=α _J1HRV ₁+ ... + α _JmHRV _m

At last according to formula L=M * HRV ^-1, can obtain L, M representes knee joint angle, and HRV representes that the user is applied to the handle retroaction vector of power on the walker, and L representes the relation between HRV and the M.

The said chaos particle swarm optimization tuning PID parameter of utilizing further is refined as:

Ratio infinitesimal analysis PID adopts ratio unit P, integral unit I and differentiation element D three parts to form, according to the error of system, through the K that sets _p, K _iAnd K _dThree parameters are controlled system:

$\mathrm{yout}\left(t\right)=K_p\mathrm{error}\left(t\right)+K_i\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}\mathrm{error}\left(j\right)+K_d[\mathrm{error}\left(t\right)-\mathrm{error}(t-1)]---\left(9\right)$

K wherein _pBe scale-up factor, K _iBe integral coefficient, K _dBe differential coefficient, error is the deviation of preset output with actual output, and u (t) is the output of PID, is again the input of controlled system simultaneously, by PID output formula $\mathrm{Yout}\left(t\right)=K_p\mathrm{Error}\left(t\right)+K_i\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}\mathrm{Error}\left(j\right)+K_d[\mathrm{Error}\left(t\right)-\mathrm{Error}(t-1)]$ Can obtain

$u(t-1)=K_p\mathrm{error}(t-1)+K_i\underset{j=0}{\overset{t-1}{\mathrm{\&Sigma;}}}\mathrm{error}\left(j\right)+K_d[\mathrm{error}(t-1)-\mathrm{error}(t-2)]---\left(10\right)$

According to:

Δu(t)＝u(t)-u(t-1)

＝K _p(error(t)-error(t-1))+K _ierror(t)+K _d(error(t)-2error(t-1)+error(t-2))

…………………………………………………………… (11)

Have:

u(t)＝Δu(t)+u(t-1)＝

u(t-1)+K _p(error(t)-error(t-1))+K _ierror(t)+K _d(error(t)-2error(t-1)+error(t-2))

………………(12)

Adopt the chaos particle swarm optimization to carry out the adaptive optimization of ratio infinitesimal analysis pid control parameter; Selecting to receive the rope space is 3 dimensions, promptly is respectively three parameters of PID controller, chooses population size m=20; The initial velocity of colony and position produce in certain spatial dimension at random, are expressed as respectively: v _i=(v _I1, v _I2, v _I3), x _i=(x _I1, x _I2, x _I3), remember that it is p that i particle searches optimal location so far _i=(p _I1, p _I2, p _I3), whole population searches to such an extent that optimal location is p up to now _Gi=(p _Gi1, p _Gi2, p _Gi3), wherein, i=1,2 ..., 20, particle swarm optimization algorithm adopts following formula that population is operated,

v _id←v _id+c ₁r ₁(p _id-x _id)+c ₂r ₂(p _gid-x _id)+c ₃r ₃(q _id-x _id) (13)

x _id←x _id+v _id (14)

Wherein, i=1,2 ..., 20; Study factor c ₁, c ₂And c ₃Be nonnegative number, general value is 0.5; r ₁, r ₂And r ₃Be the random number between [0,1], q _IdIt is the picked at random particle position;

The concrete steps that realize are:

1, confirms parameter: study factor c ₁, c ₂And c ₃And the scale N of colony, evolution number of times and chaos optimizing number of times;

2, producing N particle at random operates;

3, by formula operate particle (13) and (14);

4, to optimal location p _Gi=(p _Gi1, p _Gi2, p _Gi3) carry out chaos optimization, with p _Gid(i=1,2 ..., 20), be mapped to

Logistic equation z _I+1=μ z _i(1-z _i) i=0,1,2 ... Field of definition [0,1]; I=1,2 ..., 20, then, carry out iteration with the Logistic equation and produce Chaos Variable

M=1,2 ... Again the Chaos Variable sequence that produces

(m=1,2 ...) through inverse mapping

(m=1,2 ...) turn back to former solution space, get

(m=1,2 ...)

In former solution space each feasible solution to the Chaos Variable experience

(m=1,2 ...) calculate its adaptive value, the feasible solution p that retention property is best ^*

5, a particle of from current colony, selecting is at random used p ^*Replace;

6, if reach maximum algebraically or obtain satisfactory solution, then optimizing process finishes, otherwise returns step 3.

Characteristics of the present invention are: utilize the HRV variation prediction knee joint angle of walk helper to change; Optimize scale-up factor, differential coefficient and the integral coefficient of PID then through the chaos particle cluster algorithm; Then control the current impulse intensity of FES system, improved FES system accuracy and stability effectively.

Description of drawings

Fig. 1 handle retroaction vector (HRV) definition synoptic diagram.

Fig. 2 is based on the FES system architecture diagram of HRV.

The structured flowchart of Fig. 3 chaos particle cluster algorithm tuning PID parameter control method.

Manikin in Fig. 4 walk-aiding functional electric stimulation.

Fig. 5 experiment scene.

The result is followed the trail of in the PID control that Fig. 6 chaos particle cluster algorithm is adjusted.

Input pass angle and the actual relative error of exporting are preset in the control of Fig. 7 chaos particle cluster algorithm tuning PID parameter down.

Embodiment

Structure based on the control of the precision in the functional electrostimulation walk help of HRV The application of new technique is as shown in Figure 2; Its workflow is: at first, utilize the HRV forecasting knee joint angle of walk help process, secondly; Utilize chaos particle swarm optimization tuning PID parameter, real-time monitoring FES levels of current intensity.Its structural representation of adjusting is as shown in Figure 3, for: at first according to three decision variable K of PID _p, K _iAnd K _dThe bound of span; Confirm parameters such as population population size, search volume dimension, and the speed and the position of initialization particle colony, the fitness value that calculates each particle in the population through the corresponding relation of actual joint angles and muscle model output joint angles as appropriate evaluation function utilized then; And its fitness and optimum position fitness value itself made comparisons; And with it as the particle typical value, then in other parameters such as adjustment particle's velocity etc., change the optimum position of particle; Till stable, calculate the K that final best position promptly gets PID _p, K _iAnd K _dThree coefficients.Computing system output yout under the new PID coefficient and with the deviation of muscle model after get into the self study and the weighting coefficient self-adjusting of next step neural network again.This process finally realizes the self-adaptation on-line tuning of pid control parameter repeatedly, and is used for the FES system.

One, HRV forecasting knee joint angle model

In the walk help process, when the user under the functional electrostimulation effect, when lifting leg and taking a step; In order to support body steadiness, user's applied force on walker is then different, because varying in size of joint can make the gravity center of human body be in diverse location; It is also different then to overcome the gravity applied force; The residing planimetric position of human body also changes to some extent simultaneously, and applied force also changes to some extent for the position is tumbled then in the plane, therefore; Joint angles and user have certain relation to the walker applied force, and be as shown in Figure 4.

M＝L·HRV+wPW (1)

Wherein, M representes knee joint angle, and HRV representes that the user is applied to the handle retroaction vector of power on the walker, and L representes the relation between HRV and the M, and w representes coefficient, and W representes the center of gravity of upper arm, trunk and lower limb, and P representes the relation between three centers of gravity and the M.

In the reality, because the effect of walker, the gravity center of human body moves less, and knee joint angle then can be expressed as

M＝L·HRV (2)

Wherein, M representes knee joint angle, and HRV representes that the user is applied to the handle retroaction vector of power on the walker, and L representes the relation between HRV and the M.Shown in formula 2, confirm that L just can utilize HRV to take out the knee joint angle in the corresponding moment.

L＝M□HRV ^-1 (3)

When this patent is found the solution L, adopted the method for PLS.

Be provided with m HRV variable HRV1 ..., HRVm, p M variable, M1 ..., Mp, common i (i=1 ..., the n) data set of individual observed reading.T, U are respectively the composition that from HRV variable and M variable, extracts, and the composition that extracts here is commonly referred to the offset minimum binary factor.

Concentrate the linear combination of extracting first couple of composition T1, U1 to be from original variable:

T ₁＝ω ₁₁HRV ₁+…+ω _1mHRV _m＝ω ₁′?HRV (4)

U ₁＝v ₁₁M ₁+…+v _1pM _p＝v ₁′M (5)

ω wherein ₁=(ω ₁₁..., ω _1m) ' be model effect weight, v ₁=(v ₁₁..., v _1p) be M variable weight.For guaranteeing that T1, U1 extract the variation information of place set of variables separately as much as possible; Guarantee that simultaneously degree of correlation between the two reaches maximum; According to the character that the covariance of composition can be calculated by the inner product of the score vector of corresponding composition, the requirement of said extracted first composition is converted into asks conditional extremum to ask.

Wherein t1, u1 are the score vector of first pair of composition of being tried to achieve by sample, and HRV0, M0 are initializaing variable.Utilize method of Lagrange multipliers, the problems referred to above are converted into asks vector of unit length ω ₁And v ₁, make θ ₁=ω ₁HRV ₀' M ₀v ₁Maximum is promptly asked matrix H RV ₀' M ₀M ₀' HRV ₀Eigenwert and proper vector, its eigenvalue of maximum is θ ₁ ², corresponding unit character vector is exactly the ω that separates that is asked ₁, and v ₁By formula

Obtain.

Next sets up the equation of initializaing variable to T1

$\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="0"\; label="HRV"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">HRV</figure-callout>}}_0={\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1{\mathrm{\&alpha;}}_1^{\&prime;}+E_1\\ {\mathrm{<figure-callout\; id="0"\; label="M"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">M</figure-callout>}}_0={\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1{\mathrm{\&beta;}}_1^{\&prime;}+F_1\end{array}\right.---\left(7\right)$

Wherein the t1 meaning is the same, α ₁'=(α ₁₁..., α _1m), β ₁'=(β ₁₁..., β _1p) be the parameter vector when only a M measures t1, E1, F1 are respectively n * m and n * p residual error battle array.Can try to achieve coefficient vector α according to common least square method ₁And β ₁, α wherein ₁Become model effect load capacity.

Can not reach the precision of regression model like first composition that extracts, utilization residual error battle array E1, F1 replace X0, Y0, repeat to extract composition, and the like.Suppose finally to have extracted r composition, HRV0, M0 to the regression equation of r composition are:

$\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="0"\; label="HRV"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">HRV</figure-callout>}}_0={\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1{\mathrm{\&alpha;}}_1^{\&prime;}+...+t_r{\mathrm{\&alpha;}}_r^{\&prime;}+{\mathrm{<figure-callout\; id="0"\; label="E\; r\; M"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">E</figure-callout>}}_r& \\ M_{}{}_0={\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout>}}_1{\mathrm{\&beta;}}_1^{\&prime;}+...+t_r{{\mathrm{\&beta;}}_r}^{\&prime;}+F_r\end{array}\right.---\left(8\right)$

The first step analyze extract in the gained HRV amount composition Tk (k=1 ..., r) regression equation that the M amount is set up r composition, i.e. t are brought in linear combination into _r=ω _K1HRV ₁+ ... + ω _KmHRV _mSubstitution M _j=t ₁β _1j+ ... + t _rβ _Rj(j=1 ..., p), promptly get the regression equation M of standardized variable _j=α _J1HRV ₁+ ... + α _JmHRV _m

Based on formula 3, can obtain L at last.

Two, chaos particle cluster algorithm tuning PID parameter control

PID is made up of ratio unit P, integral unit I and differentiation element D three parts, according to the error of system, through the K that sets _p, K _iAnd K _dThree parameters are controlled system.

$\mathrm{yout}\left(t\right)=K_p\mathrm{error}\left(t\right)+K_i\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}\mathrm{error}\left(j\right)+K_d[\mathrm{error}\left(t\right)-\mathrm{error}(t-1)]---\left(9\right)$

K wherein _pBe scale-up factor, K _iBe integral coefficient, K _dBe differential coefficient, error is the deviation of preset output with actual output, and u (t) is the output of PID, is again the input of controlled system simultaneously.

Can obtain by PID output formula (1)

$u(t-1)=K_p\mathrm{error}(t-1)+K_i\underset{j=0}{\overset{t-1}{\mathrm{\&Sigma;}}}\mathrm{error}\left(j\right)+K_d[\mathrm{error}(t-1)-\mathrm{error}(t-2)]---\left(10\right)$

According to:

Δu(t)＝u(t)-u(t-1)

＝K _p(error(t)-error(t-1))+K _ierror(t)+K _d(error(t)-2error(t-1)+error(t-2))

……………………………………………………………(11)

Have:

u(t)＝Δu(t)+u(t-1)＝

u(t-1)+K _p(error(t)-error(t-1))+K _ierror(t)+K _d(error(t)-2error(t-1)+error(t-2))

………………(12)

The present invention adopts chaos particle swarm algorithm to carry out the adaptive optimization of pid control parameter; Selecting to receive the rope space is 3 dimensions, promptly is respectively three parameters of PID controller, chooses population size m=20; The initial velocity of colony and position produce in certain spatial dimension at random, are expressed as respectively: v _i=(v _I1, v _I2, v _I3), x _i=(x _I1, x _I2, x _I3), remember that it is p that i particle searches optimal location so far _i=(p _I1, p _I2, p _I3), whole population searches to such an extent that optimal location is p up to now _Gi=(p _Gi1, p _Gi2, p _Gi3).Wherein, i=1,2 ..., 20, particle swarm optimization algorithm adopts following formula that population is operated.

v _id←v _id+c ₁r ₁(p _id-x _id)+c ₂r ₂(p _gid-x _id)+c ₃r ₃(q _id-x _id) (13)

x _id←x _id+v _id (14)

Wherein, i=1,2 ..., 20; Study factor c ₁, c ₂And c ₃Be nonnegative number, general value is 0.5; r ₁, r ₂And r ₃Be the random number between [0,1], q _IdIt is the picked at random particle position.

The concrete steps of its realization are:

1, confirms parameter: study factor c ₁, c ₂And c ₃And the scale N of colony, evolution number of times and chaos optimizing number of times;

2, producing N particle at random operates;

3, by formula operate particle (13) and (14);

4, to optimal location p _Gi=(p _Gi1, p _Gi2, p _Gi3) carry out chaos optimization, with p _Gid(i=1,2 ..., 20), be mapped to Logistic equation z _I+1=μ z _i(1-z _i) i=0,1,2 ... Field of definition [0,1];

I=1,2 ..., 20, then, carry out iteration with the Logistic equation and produce Chaos Variable M=1,2 ..., again the Chaos Variable sequence that produces (m=1,2 ...) through inverse mapping

(m=1,2 ...) turn back to former solution space, get

(m=1,2 ...)

In former solution space each feasible solution to the Chaos Variable experience

(m=1,2 ...) calculate its adaptive value, the feasible solution p that retention property is best ^*

5, a particle of from current colony, selecting is at random used p ^*Replace.

6, if reach maximum algebraically or obtain satisfactory solution, then optimizing process finishes, otherwise returns step 3.

Three, experimental program

Experimental provision adopts the walker system of wireless transmission and the Parastep functional electric stimulation system that U.S. SIGMEDICS company produces, and this system comprises microprocessor and boost pulse generation circuit, contains six stimulation channels, powered battery.Experiment content is: utilize the FES system that the relevant muscle group of lower limb is stimulated; Make the experimenter according to predetermined actions, record is applied to HRV on the walker at first through being installed in voltage signal and the knee joint angle movement locus that foil gauge (BX3506AA) network of electrical bridge changes into that lead of 12 on the walker simultaneously.Require the experimenter healthy, no lower limb muscles, bone illness, impassivity illness and severe cardiac pulmonary disease.The experimenter sits idly before walker during experiment, and stimulating electrode is fixed in corresponding position, and when not applying electro photoluminescence, it is light that the experimenter keeps.The FES experiment scene is as shown in Figure 5.The electric stimulation pulse sequence adopts classical Lilly waveform, and pulsed frequency is 25Hz, pulsewidth 150 μ s, and pulse current is adjustable in 0～120m scope.In the experiment, write down HRV in real time and can adjust stimulus intensity to change the knee joint angle that produces by stimulating through changing the pulse current size.Before the experiment, set the knee joint angle movement locus of expectation, utilize the measurement of angle meter to detect the knee joint subtended angle in real time in the experiment and change.The experimental data sampling rate is 128Hz, and the data recording duration is 60s.

Beneficial effect

The new algorithm of chaos particle swarm algorithm tuning PID parameter is calculated the FES pulse current amplitude and is adjusted, the knee joint angle that the FES effect is produced move the movement locus of expection.The result is followed the trail of in the PID control that Fig. 6 adjusts for chaos particle swarm algorithm.Red line representes that desired movement track, blue line are actual output joint angles among the figure.The X axle is the time, and the Y axle is the motion of knee joint angle.For more clearly observing the departure of chaos particle cluster algorithm tuning PID; Shown in the relative error of preset input knee joint angle and actual knee joint angle under Fig. 7 chaos particle swarm Tuning PID Controller; Can find out that then error all within 5%, can reach accurate control.

Purport of the present invention is the precision control method that proposes a kind of new FES; Utilize the error of knee joint angle and the joint angles of actual knee joint angle prediction of the HRV parameter prediction of walker; Through scale-up factor, integral coefficient and the differential coefficient of chaos swarm optimization algorithm PID, the accurately stable then systematically strength of current of FES of controlling in real time.This invention can improve FES system accuracy and stability effectively, and obtains considerable social benefit and economic benefit.Optimum implementation is intended and is adopted patent transfer, technological cooperation or product development.

Claims (2)

1. a functional electrostimulation pid parameter double source characteristic fusion particle swarm setting method is characterized in that, comprises the following steps:

At first, utilize the handle retroaction vector HRV forecasting knee joint angle of walk help process;

Secondly, utilize the chaos particle swarm optimization ratio infinitesimal analysis pid parameter of adjusting, real-time monitoring FES levels of current intensity, the flow process of adjusting is: at first according to three decision variable K of ratio infinitesimal analysis pid parameter _p, K _iAnd K _dThe bound of span; Confirm to comprise the parameter of population population size, search volume dimension; And the speed and the position of initialization particle colony; Utilize then through actual joint angles and muscle model and export the fitness value of the corresponding relation of joint angles as each particle in the appropriate evaluation function calculating population, corresponding relation is L=M * HRV ^-1, wherein, M representes knee joint angle; HRV representes that the user is applied to the handle retroaction vector of power on the walker, and L representes the relation between HRV and the M, adopts the method for PLS to confirm muscle model output joint angles; And the fitness value and the optimum position fitness value of this particle made comparisons; And as the particle typical value, and then adjustment particle's velocity and other parameters change the optimum position of particle with the optimum position of the fitness value of this particle and this particle itself; Till stable, calculate three decision variable K that final best position promptly gets ratio infinitesimal analysis pid parameter _p, K _iAnd K _d, computing function property electro photoluminescence FES system output yout under the new ratio infinitesimal analysis pid parameter and with the deviation of muscle model output joint angles after get into the self study of next step chaos particle swarm optimization neural network and three decision variable K of pid parameter at that time again _p, K _iAnd K _dThe self-adjusting of weighting coefficient, this process finally realizes the self-adaptation on-line tuning of ratio infinitesimal analysis pid parameter repeatedly, and is used for functional electrostimulation FES system.

2. a kind of functional electrostimulation pid parameter double source characteristic fusion particle swarm setting method according to claim 1 is characterized in that, muscle model output joint angles is the method that adopts PLS, that is:

Be provided with m HRV variable HRV1 ..., HRVm, p M variable, M1 ..., Mp, the data set of common i observed reading, i=1 ..., n; T, U are respectively the composition that from HRV variable and M variable, extracts, and are called the offset minimum binary factor, concentrate the linear combination of extracting first couple of composition T1, U1 to be from original variable:

T ₁＝ω ₁₁HRV ₁+…+ω _1mHRV _m＝ω′ ₁HRV (4)

U ₁＝v ₁₁M ₁+…+v _1pM _p＝v′ ₁M (5)

ω wherein ₁=(ω ₁₁..., ω _1m) ' be model effect weight, v ₁=(v ₁₁..., v _1p) ' be M variable weight is converted into the requirement of first pair of composition of said extracted and asks constrained extremal problem:

Wherein t1, u1 are the score vector of first pair of composition of being tried to achieve by sample, and HRV0, M0 are initializaing variable, utilize method of Lagrange multipliers, and the problems referred to above are converted into asks vector of unit length ω ₁And v ₁, make θ ₁=ω ' ₁HRV ' ₀M ' ₀v ₁' maximum is promptly asked matrix H RV ' ₀M ₀M ' ₀HRV ₀Eigenwert and proper vector, its eigenvalue of maximum is θ ₁ ², corresponding unit character vector is exactly the ω that separates that is asked ₁, and v ₁By formula Obtain;

Next sets up the equation of initializaing variable to T1

Wherein the t1 meaning is the same, α ' ₁=(α ₁₁..., α _1m), β ' ₁=(β ₁₁..., β _1p) be the parameter vector during a M variable t1 only, E1, F1 are respectively n * m and n * p residual error battle array, try to achieve coefficient vector α according to common least square method ₁And β ₁, α wherein ₁Be called model effect load capacity;

Can not reach the precision of regression model like the first pair of composition that extracts, utilization residual error battle array E1, F1 replace HRV ₀, M ₀, repeat to extract composition, and the like, supposing finally to have extracted r composition, HRV0, M0 to the regression equation of r composition are:

The first step analyze extract in the gained HRV amount composition Tk (k=1 ..., r) regression equation that the M variable is set up r composition, i.e. t are brought in linear combination into _r=ω _K1HRV ₁+ ... + ω _KmHRV _mSubstitution M _j=t ₁β _1j+ ... T _rβ _Rj(j=1 ..., p), promptly get the regression equation M of standardized variable _j=α _J1HRV ₁+ ... + α _JmHRV _m

At last according to formula L=M * HRV ^-1, can obtain L.

3 .A kind of functional electrostimulation pid parameter double source characteristic fusion particle swarm setting method according to claim 1 is characterized in that, utilizes chaos particle swarm optimization tuning PID parameter, further is refined as:

Ratio infinitesimal analysis PID adopts ratio unit P, integral unit I and differentiation element D three parts to form, according to the error of functional electric stimulation system, through the K that sets _p, K _iAnd K _dThree parameters are controlled functional electric stimulation system:

K wherein _pBe scale-up factor, K _iBe integral coefficient, K _dBe differential coefficient, error is the deviation of preset output with actual output, and u (t) is the output of PID, is again the input of controlled system simultaneously, by PID output formula

can obtain

According to:

Δu(t)＝u(t)-u(t-1)

＝K _p(error(t)-error(t-1))+K _ierror(t)+K _d(error(t)-2error(t-1)+error(t-2))……………………………………………………………(11)

Have:

u(t)＝Δu(t)+u(t-1)＝

u(t-1)+K _p(error(t)-error(t-1))+K _ierror(t)+K _d(error(t)-2error(t-1)+error(t-2))

………………(12)

Adopt the chaos particle swarm optimization to carry out the adaptive optimization of ratio infinitesimal analysis pid control parameter; Selecting the search volume is 3 dimensions, promptly is respectively three parameters of PID controller, chooses population size m=20; The initial velocity of colony and position produce in certain spatial dimension at random, are expressed as respectively: v _i=(v _I1, v _I2, v _I3), x _i=(x _I1, x _I2, x _I3), remember that it is p that i particle searches optimal location so far _i=(p _I1, p _I2, p _I3), whole population searches to such an extent that optimal location is p up to now _Gi=(p _Gi1, p _Gi2, p _Gi3), wherein, i=1,2 ..., 20, particle swarm optimization algorithm adopts following formula population to be operated v _Id← v _Id+ c ₁r ₁(p _Id-x _Id)+c ₂r ₂(p _Gid-x _Id)+c ₃r ₃(q _Id-x _Id) (13)

x _id←x _id+v _id (14)

Wherein, study factor c ₁, c ₂And c ₃Be nonnegative number, value is 0.5; r ₁, r ₂And r ₃Be the random number between [0,1], q _IdBe the picked at random particle position, d representes the coding method of ion, d=10;

The concrete steps that realize are:

1, confirms parameter: study factor c ₁, c ₂And c ₃And the scale N of colony, evolution number of times and chaos optimizing number of times;

2, producing N particle at random operates;

3, by formula operate particle (13) and (14);

4, to optimal location p _Gi=(p _Gi1, p _Gi2, p _Gi3) carry out chaos optimization, with p _GidBe mapped to Logistic equation z _I+1=μ z _i(1-z _i), z _I+1Be the value of Chaos Variable, its field of definition is [0,1], and μ is the Logistic controlled variable, the control chaos state;

Then, carry out iteration with the Logistic equation and produce Chaos Variable

M=1,2 ..., again the Chaos Variable sequence that produces

Through inverse mapping

Turn back to former solution space,

a _IdEach minimum parameter in the expression pid parameter, b _IdEach maximum parameter in the expression pid parameter;

In former solution space each feasible solution to the Chaos Variable experience

Calculate its adaptive value, the feasible solution p that retention property is best ^*

5, a particle of from current colony, selecting is at random used p ^*Replace;

6, if reach maximum iteration time or obtain satisfactory solution, then optimizing process finishes, otherwise returns step 3.

Images