CN101837164A - Double source feature fusion ant colony tuning method for PID (Proportion Integration Differention) parameter in functional electro-stimulation - Google Patents

Abstract

The invention relates to the field of equipment for limb rehabilitation by electric pulse stimulation. In order to realize accurate, stable and real-time control of the current intensity of the FES system and effectively improve the accuracy and stability of the FES system, the technical solution adopted in the present invention is: a dual-source feature fusion ant colony tuning method for PID parameters in functional electrical stimulation, including The following steps: first, use the muscle model HRV of the walking aid process to predict the knee joint angle; second, use the ant colony algorithm to adjust the PID parameters, and adjust the FES current level in real time; use the ant colony algorithm to control the PID parameters, if the expected goal is not reached Continue to optimize; calculate the system output and its deviation from the muscle model HRV under the new PID coefficient, and then enter the next step of ant colony algorithm self-learning and weight coefficient self-adjustment; repeat the previous step to realize the self-adaptation of PID control parameters Online setting, and used in FES system. The invention is mainly applied to setting PID parameters.

Description

The double source feature fusion ant colony setting method of pid parameter in the functional electric stimulation

Technical field

The present invention relates to carry out the instrument field of limb rehabilitating, especially the double source feature fusion ant colony setting method of pid parameter in the functional electric stimulation with electric pulse stimulation.

Background technology

(Functional Electrical Stimulation is to stimulate limb motion muscle group and peripheral nervous thereof by current pulse sequence FES) to functional electric stimulation, 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, at 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 electricity irritation peroneal nerve to correct the gait of hemiplegic patient's drop foot first, has started the new way that functional electric stimulation 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 finished the key problem in technology that the intended function action is still FES with assurance electricity irritation 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, automatically adjust FES stimulus intensity and time sequence parameter with closed loop control, thereby improved the accuracy and the stability of FES system greatly, but now effective control method is still among exploring.

Handle retroaction vector (handle reactions vector, HRV) be according in the process of standing and walking under the walker help, in fact the effectiveness that walker offers the patient can be divided into clear and definite independently 3 parts: 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 concentrfated load to the effect of walker is synthetic in the walking process of standing, represent with two mechanics vectors at handle mid point cross section centre of form place respectively, as shown in Figure 1, vector is at x, y, durection component on z axle size with joint efforts can characterize the patient respectively by trying hard to recommend into that walker obtained, dynamic balance and power support level.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 electricity irritation 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 electricity irritation action effect can accurately be finished the intended function action, and prevents muscle fatigue.

Ratio calculus (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 complexity and the time variation operating environment of muscle, because good stability, the reliable operation of PID have still obtained in the functional electric stimulation field using widely at present.The PID core technology is accurate determine 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 deficiencies in the prior art, the double source feature fusion ant colony setting method of pid parameter in a kind of functional electric stimulation is provided, can accurately stablize and control systematically current intensity 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 solution used in the present invention is: the double source feature fusion ant colony setting method of pid parameter in the functional electric stimulation comprises the following steps:

At first, utilize the muscle model HRV forecasting knee joint angle of walk help process;

Secondly, utilize the ant group algorithm pid parameter of adjusting, real-time monitoring FES levels of current intensity, its flow process of adjusting is: head

Elder generation is according to three decision variable K of PID _p, K _iAnd K _dThe bound of span, determine to comprise the parameter of ant group population size, search volume dimension, and it is encoded, utilize then by actual joint angles and export the fitness value of the corresponding relation of joint angles as appropriate evaluation function calculating with muscle model HRV; Adopt ant group algorithm that pid parameter is controlled, promptly determine the parameter setting of ant group algorithm, utilize the Formica fusca random search to make its variable optimize the K of PID _p, K _iAnd K _dThree coefficients utilize fitness function to regulate at every turn explore path and judge whether to reach goal-selling of Formica fusca, as reach goal-selling, calculate the K that final best position promptly gets PID _p, K _iAnd K _dThree coefficients, as do not reach re-set target continuation optimizing, up to reaching goal-selling; Computing system output yout under the new PID coefficient and with the deviation of muscle model HRV after enter next step ant group algorithm again self study and weight coefficient self-adjusting;

Previous step finally realizes the self adaptation on-line tuning of pid control parameter repeatedly, and is used for the FES system.

The corresponding relation of described actual joint angles and muscle model HRV output joint angles is:

L＝M×HRV ^-1 (3)

M represents knee joint angle, and HRV represents that user is applied to the handle retroaction vector of power on the walker, and L represents the relation between HRV and the M, adopts the method for PLS to find the solution L:

Be provided with m HRV variable HRV1 ..., HRVm, p M variable, M1 ..., Mp, i (i=1 altogether,, the n) data set of individual observation, T, U are respectively the composition that extracts from HRV variable and M variable, 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 unit vector ω ₁And v ₁, make θ ₁=ω ₁' HRV ₀' M ₀v ₁Maximum is promptly asked matrix H RV ₀' M ₀M ₀' HRV ₀Eigenvalue and characteristic 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="HDA0000021781220000031.png"\; state="\{\{state\}\}">HRV</figure-callout>}}_0=t_1{\mathrm{\&alpha;}}_1^{\&prime;}+E_1\\ {\mathrm{<figure-callout\; id="0"\; label="M"\; filenames="HDA0000021781220000031.png"\; state="\{\{state\}\}">M</figure-callout>}}_0=t_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 method of least square ₁And β ₁, α wherein ₁Become model effect load capacity;

Can not reach the precision of regression model as 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="HDA0000021781220000031.png"\; state="\{\{state\}\}">HRV</figure-callout>}}_0=t_1{\mathrm{\&alpha;}}_1^{\&prime;}+...+t_r{\mathrm{\&alpha;}}_r^{\&prime;}+{\mathrm{<figure-callout\; id="0"\; label="E\; r\; M"\; filenames="HDA0000021781220000031.png"\; state="\{\{state\}\}">E</figure-callout>}}_r& \\ M_{}{}_0=t_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

According to formula (3), can obtain L at last.

The described coding is that according to knee joint angle and current-mode, and three parameters of actual error situation setting PID are 5 position effective digitals, wherein K _pPreceding 2 of arithmetic point, behind the arithmetic point 3; K _iAnd K _dPreceding 1 of arithmetic point, behind the arithmetic point 4, specific coding such as following formula:

K _p＝y _1，j×10 ¹+y _2，j×10 ⁰+y _3，j×10 ^-1+y _4，j×10 ^-2+y _5，j×10 ^-3

K _i＝y _6，j×10 ⁰+y _7，j×10 ^-1+y _8，j×10 ^-2+y _9，j×10 ^-3+y _10，j×10 ^-4

K _d＝y _6，j×10 ⁰+y _7，j×10 ^-1+y _8，j×10 ^-2+y _9，j×10 ^-3+y _10，j×10 ^-4

And the concrete function of overshoot equal error design adaptation function according to output is:

fit＝α1×σ+β1×t+c×error (9)

Wherein, σ is an overshoot, and t is the rise time, and error is the relative error of output joint angles and default joint angles, α 1=0.1, β 1=0.8, c=2;

The concrete workflow of ant group algorithm Tuning PID Controller device parameter is:

Step1: parameter initialization makes time t=0 and cycle-index N _Max=0, maximum cycle N is set _Cmax, m Formica fusca placed starting point;

Step2: Formica fusca number and cycle-index are set;

Step3: Formica fusca random search, after the end of once creeping, determine the actual input variable of the selected conduct of which characteristic variable, revise the taboo list index, after promptly choosing Formica fusca is moved to new element, and this element is moved in the taboo table of this Formica fusca individuality;

Step4: calculate the probability that the Formica fusca individuality calculates according to the state transition probability formula, select element according to this probability;

Step5:, change Step3, otherwise be Step6 if the Formica fusca element has not traveled through;

Step6: the pheromone concentration that the plain concentration of lastest imformation is divided the high characteristic variable of accuracy is enhanced, and next time can be selected with bigger probability when searching for;

Step7: satisfy and finish to regulate the end of adjusting.

Characteristics of the present invention are: utilize the HRV variation prediction knee joint angle of walking aid to change, pass through proportionality coefficient, differential coefficient and the integral coefficient of ant group group algorithm optimization PID then, 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 sketch map.

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

Fig. 3 ant group group algorithm structured flowchart of pid parameter control method of adjusting.

Anthropometric dummy in Fig. 4 walk-aiding functional electric stimulation.

Fig. 5 ant group algorithm pid parameter coding sketch map of adjusting.

Fig. 6 experiment scene.

The result is followed the trail of in the PID control that Fig. 7 ant group group algorithm is adjusted.

The specific embodiment

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, by proportionality coefficient, integral coefficient and the differential coefficient of ant group algorithm optimization PID, the accurately stable then systematically current intensity 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.

Based on the structure of the precision in the functional electric stimulation walk help of HRV control The application of new technique as shown in Figure 2, its workflow is: at first, utilize the HRV forecasting knee joint angle of walk help process, secondly, utilize the ant group algorithm pid parameter of adjusting, real-time monitoring FES levels of current intensity.It adjusts structural representation as shown in Figure 3, for: at first according to three decision variable K of PID _p, K _iAnd K _dThe bound of span, determine parameters such as ant group population size, search volume dimension, and it is encoded, utilize the fitness value that calculates as appropriate evaluation function by the corresponding relation of actual joint angles and muscle model output joint angles then, and the parameter setting of definite ant group algorithm, utilize Formica fusca to receive rope at random and make its variable optimize the K of PID _p, K _iAnd K _dThree coefficients utilize fitness function to regulate at every turn explore path and judge whether to reach goal-selling of Formica fusca.As reach goal-selling, calculate the K that final best position promptly gets PID _p, K _iAnd K _dThree coefficients, as do not reach re-set target continuation optimizing, up to reaching goal-selling.Computing system output yout under the new PID coefficient and with the deviation of muscle model after enter next step ant group algorithm again self study and weight coefficient self-adjusting.This process finally realizes the self adaptation on-line tuning of pid control parameter repeatedly, and is used for the FES system.

1HRV forecasting knee joint angle model

In the walk help process, when user under the functional electric stimulation effect, when lifting lower limb and taking a step, in order to support body steadiness, user 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 plan-position of human body also changes to some extent simultaneously, 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, as shown in Figure 4.

M＝L·HRV+wPW (1)

Wherein, M represents knee joint angle, and HRV represents that user is applied to the handle retroaction vector of power on the walker, and L represents the relation between HRV and the M, and w represents coefficient, and W represents the center of gravity of upper arm, trunk and lower limb, and P represents 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 represents knee joint angle, and HRV represents that user is applied to the handle retroaction vector of power on the walker, and L represents the relation between HRV and the M.Shown in formula 2, determine that L just can utilize HRV to take out the knee joint angle in the corresponding moment.

L＝M×HRV ^-1 (3)

When the present invention finds 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 observation.T, U are respectively the composition that extracts from HRV variable and M variable, the composition of Ti Quing is commonly referred to the offset minimum binary factor here.

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 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, the problems referred to above are converted into asks unit vector ω ₁And v ₁, make θ ₁=ω ₁' HRV ₀' M ₀v ₁Maximum is promptly asked matrix H RV ₀' M ₀M ₀' HRV ₀Eigenvalue and characteristic 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="HDA0000021781220000031.png"\; state="\{\{state\}\}">HRV</figure-callout>}}_0=t_1{\mathrm{\&alpha;}}_1^{\&prime;}+E_1\\ {\mathrm{<figure-callout\; id="0"\; label="M"\; filenames="HDA0000021781220000031.png"\; state="\{\{state\}\}">M</figure-callout>}}_0=t_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 method of least square ₁And β ₁, α wherein ₁Become model effect load capacity.

Can not reach the precision of regression model as 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="HDA0000021781220000031.png"\; state="\{\{state\}\}">HRV</figure-callout>}}_0=t_1{\mathrm{\&alpha;}}_1^{\&prime;}+...+t_r{\mathrm{\&alpha;}}_r^{\&prime;}+{\mathrm{<figure-callout\; id="0"\; label="E\; r\; M"\; filenames="HDA0000021781220000031.png"\; state="\{\{state\}\}">E</figure-callout>}}_r& \\ M_{}{}_0=t_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

According to formula 3, can obtain L at last.

Ant group algorithm is at first encoded to three parameters of PID to the control of pid parameter, and according to knee joint angle and current-mode, and three parameters that situation such as actual error is set PID are 5 position effective digitals, wherein K _pPreceding 2 of arithmetic point, behind the arithmetic point 3; K _iAnd K _dPreceding 1 of arithmetic point, behind the arithmetic point 4.Its specific coding sketch map as shown in Figure 5.

As following formula

K _p＝y _1，j×10 ¹+y _2，j×10 ⁰+y _3，j×10 ^-1+y _4，j×10 ^-2+y _5，j×10 ^-3

K _i＝y _6，j×10 ⁰+y _7，j×10 ^-1+y _8，j×10 ^-2+y _9，j×10 ^-3+y _10，j×10 ^-4

K _d＝y _6，j×10 ⁰+y _7，j×10 ^-1+y _8，j×10 ^-2+y _9，j×10 ^-3+y _10，j×10 ^-4

And the concrete function of overshoot equal error design adaptation function according to output is:

fit＝α1×σ+β1×t+c×error (9)

Wherein, σ is an overshoot, and t is the rise time, and error is the relative error of output joint angles and default joint angles, α 1=0.1, β 1=0.8, c=2.

The concrete workflow of ant group algorithm Tuning PID Controller device parameter is:

Step1: parameter initialization.Make time t=0 and cycle-index N _Max=0, maximum cycle N is set _Cmax, m Formica fusca placed starting point.

Step2: Formica fusca number and cycle-index are set

Step3: the Formica fusca random search, after the end of once creeping, determine the actual input variable of the selected conduct of which characteristic variable, revise the taboo list index, after promptly choosing Formica fusca is moved to new element, and this element is moved in the taboo table of this Formica fusca individuality

Step4: calculate the probability that the Formica fusca individuality calculates according to the state transition probability formula, select element according to this probability

Step5:, change Step3, otherwise be Step6 if the Formica fusca element has not traveled through

Step6: the pheromone concentration that the plain concentration of lastest imformation is divided the high characteristic variable of accuracy is enhanced, and next time can be selected with bigger probability when searching for

Step7: satisfy and finish to regulate the end of adjusting.

The control that 2 ant group group algorithms are adjusted pid parameter

PID is made up of ratio unit P, integral unit I and differentiation element D three parts, according to the error of system, by 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 proportionality coefficient, K _iBe integral coefficient, K _dBe differential coefficient, error is the deviation of default 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 ant group algorithm to carry out the adaptive optimization of pid control parameter, as a combination, utilizes three parameters of PID the ant group algorithm optimizing to solve this combinatorial problem.Ant group algorithm is a kind of novel bionic Algorithm that comes from the Nature biological world, when finding the solution optimization problem with ant group algorithm, at first optimization problem is transformed in order to find the solution shortest route problem.Every Formica fusca is from initial contact N ₀₀, N ₀₁... N _0nSet out, N in proper order passes by ₁, N ₂..., a wherein child node, up to destination node N _K0, N _K1... N _KnForm path (N _0tN _1t... N _Kt), t ∈ [0,1,2 ... 9].A binary feasible solution can be represented in its path.Following feature is arranged during each Formica fusca visit city:

The state transformation rule: the state transformation rule that ant group algorithm uses is the rule of ratio at random that proposes based on the TSP problem, and it provides the probability that the Formica fusca k that is positioned at city i selects to move to city j,

τ wherein _Ij(i j) is (i, fitness j), η _Ij(i, j)=(10-|y (i)-y (i) ^*|)/10, during y (i) ant group hunting in the value at i place, y (i) ^*When searching for for last time in the value at i place.α is the relative significance level of residual risk, the relative significance level that β is expected value.

In ant group algorithm, selection mode is

Wherein, q is for being evenly distributed on a random number on [0,1], q ₀Be the parameter on [0,1].

Overall situation update rule: ant algorithm has different update algorithm, the overall situation that ant group system adopts is upgraded principle, only allowing the Formica fusca release pheromone of globally optimal solution, is the neighborhood that mainly concentrates on the best path of being found out till the current circulation for the search that makes Formica fusca like this.

τ _ij(i，j)←(1-ρ)□τ _ij(i，j)+ρ·Δτ _ij(i，j) (14)

Wherein ρ is that information is counted volatility coefficient, L _GbBe the global optimum path of finding so far

Local updating information: every Formica fusca is set up the renewal that the plain mark of the information of carrying out number is also arranged in the process of separating

τ _ij(i，j)←(1-γ)□τ _ij(i，j)+γ·Δτ _ij(i，j) (16)

γ ∈ [0,1] wherein.

Ant group algorithm is at first encoded to three parameters of PID to the control of pid parameter, and according to knee joint angle and current-mode, and three parameters that situation such as actual error is set PID are 5 position effective digitals, wherein K _pPreceding 2 of arithmetic point, behind the arithmetic point 3; K _iAnd K _dPreceding 1 of arithmetic point, behind the arithmetic point 4.Its specific coding sketch map as shown in Figure 5.

As following formula

K _p＝y _1，j×10 ¹+y _2，j×10 ⁰+y _3，j×10 ^-1+y _4，j×10 ^-2+y _5，j×10 ^-3

K _i＝y _6，j×10 ⁰+y _7，j×10 ^-1+y _8，j×10 ^-2+y _9，j×10 ^-3+y _10，j×10 ^-4

K _d＝y _6，j?×10 ⁰+y _7，j×10 ^-1+y _8，j?×10 ^-2+y _9，j?×10 ^-3+y _10，j×10 ^-4

The concrete workflow of ant group algorithm Tuning PID Controller device parameter is:

Step1: parameter initialization.Make time t=0 and cycle-index N _Max=0, maximum cycle N is set _Cmax, m Formica fusca placed starting point.

Step2: Formica fusca number and cycle-index are set

Step3: the Formica fusca random search, after the end of once creeping, determine the actual input variable of the selected conduct of which characteristic variable, revise the taboo list index, after promptly choosing Formica fusca is moved to new element, and this element is moved in the taboo table of this Formica fusca individuality

Step4: calculate the probability that the Formica fusca individuality calculates according to the state transition probability formula, select element according to this probability

Step5:, change Step3, otherwise be Step6 if the Formica fusca element has not traveled through

Step6: the pheromone concentration that the plain concentration of lastest imformation is divided the high characteristic variable of accuracy is enhanced, and next time can be selected with bigger probability when searching for

Step7: satisfy and finish to regulate the end of adjusting.

3 experimental programs

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, battery powered.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 by being installed in voltage signal and the knee joint angle movement locus that foil gauge (BX350-6AA) network of electrical bridge changes into that lead of 12 on the walker simultaneously.Require the experimenter healthy, no lower limb muscles, skeleton illness, impassivity illness and severe cardiac pulmonary disease.Before the experimenter sits on walker during experiment,

Stimulating electrode is fixed in corresponding position, and when not applying electricity irritation, it is light that the experimenter keeps.The FES experiment scene as shown in Figure 5.The electric stimulation pulse sequence adopts classical Lilly waveform, and pulse frequency is 25Hz, pulsewidth 150 μ s, and pulse current is adjustable in 0～120m scope.In the experiment, write down in real time HRV with

And can adjust stimulus intensity to change the knee joint angle that produces by stimulating by changing the pulse current size.Before the experiment, set the knee joint angle movement locus of expectation, utilize the angular surveying meter to detect the knee joint subtended angle in real time in the experiment and change.The experimental data sample rate is 128Hz, and the data record duration is 60s.

Beneficial effect

The adjust new algorithm of pid parameter of ant group algorithm 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.Fig. 7 follows the trail of the result for the PID control that ant group algorithm is adjusted.Red line represents that desired movement track, blue line are actual output joint angles among the figure.X-axis is the time, and Y-axis is the motion of knee joint angle.For more clearly observing the departure that ant group algorithm is adjusted PID, shown in the relative error of default input knee joint angle and actual knee joint angle under the ant group algorithm Tuning PID Controller, then error can reach accurate control all within 5% as can be seen.

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, by proportionality coefficient, integral coefficient and the differential coefficient of ant group algorithm optimization PID, the accurately stable then systematically current intensity 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 intends adopting patent transfer, technological cooperation or product development.

Claims (4)

1. a dual-source feature fusion ant colony tuning method of PID parameters in functional electrical stimulation, is characterized in that, comprises the following steps: First, the knee joint angle is predicted using the muscle model HRV of the walking aid process; Secondly, use the ant colony algorithm to tune the PID parameters, and adjust the FES current level _{in real time} . scale, search space dimension parameters, and encode them, and then use the corresponding relationship between the actual joint angle and the muscle model HRV output joint angle as the fitness value calculated by the moderate evaluation function; use the ant colony algorithm to control the PID parameters , that is to determine the parameter setting of the ant colony algorithm, use the ant random search to make its variables optimize the three coefficients K _p , K _i and K _d of the PID, use the fitness function to adjust the ants' search path each time and judge whether the preset goal is reached, If the preset goal is reached, calculate the final best position to get the three coefficients of PID K _p , K _i and K _d , if the expected goal is not reached, continue to optimize until the preset goal is reached; calculate under the new PID coefficient After the system outputs yout and its deviation from the muscle model HRV, it enters the next step of ant colony algorithm self-learning and weighting coefficient self-adjustment; Repeat the previous step, and finally realize the adaptive online tuning of PID control parameters, and use it in the FES system. 2. the dual-source feature fusion ant colony tuning method of PID parameter in a kind of functional electrical stimulation according to claim 1, is characterized in that, the corresponding relation of actual joint angle and muscle model HRV output joint angle is: L＝M×HRV ^-1 (3) M represents the knee joint angle, HRV represents the handle reaction vector of the force exerted by the user on the walker, L represents the relationship between HRV and M, and the partial least squares regression method is used to solve L: There are m HRV variables HRV1, ..., HRVm, p M variables, M1, ..., Mp, a data set with a total of i (i=1, ..., n) observed values, T and U are respectively from HRV variables and The components extracted from the M variable, the linear combination of the first pair of components T1 and U1 extracted from the original variable set is: T ₁ =ω ₁₁ HRV ₁ +...+ω _1m HRV _m =ω′ ₁ HRV (4) U ₁ =v ₁₁ M ₁ +...+v _1p M _p =v' ₁ M (5) Where ω ₁ =(ω ₁₁ ,...,ω _1m )' is the weight of the model effect, v ₁ =(v ₁₁ ,...,v _1p )' is the weight of the M variable. Value question:

Among them, t1 and u1 are the score vectors of the first pair of components obtained from the samples, HRVO and MO are the initial variables, and using the Lagrange multiplier method, the above problem is transformed into finding unit vectors ω ₁ and v ₁ , so that θ ₁ =ω′ ₁ HRV′ ₀ M ₀ v ₁ is the largest, that is to find the eigenvalues and eigenvectors of the matrix HRV′ ₀ M ₀ M′ ₀ HRV ₀ , its largest eigenvalue is θ ₁ ² , and the corresponding unit eigenvector is the desired The solution of ω ₁ , and v ₁ is given by the formula get; Secondly, establish the equation of the initial variable to T1 $\left\{\begin{array}{c}{\mathrm{<span\; class="notranslate">\; HRV</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; E.</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\\ {\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\end{array}\right.<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 7</span>}<span\; class="notranslate">\; )</span>$ Where t1 has the same meaning as before, α′ ₁ =(α ₁₁ ,…,α _1m ), β′ ₁ =(β ₁₁ ,…,β _1p ) is the parameter vector when there is only one M quantity t1, E1 and F1 are n respectively ×m and n×p residual matrix, the coefficient vectors α ₁ and β ₁ can be obtained according to the ordinary least squares method, where α ₁ becomes the model effect load; If the extracted first component cannot reach the accuracy of the regression model, use the residual matrix E1 and F1 to replace XO and YO, repeat the extraction of components, and so on, assuming that r components are finally extracted, the regression equation of HRVO and MO for r components for: $\left\{\begin{array}{c}{\mathrm{<span\; class="notranslate">\; HRV</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; r</span>}}{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{\mathrm{<span\; class="notranslate">\; r</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; E.</span>}}_{\mathrm{<span\; class="notranslate">\; r</span>}}\\ {\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; r</span>}}{{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}_{\mathrm{<span\; class="notranslate">\; r</span>}}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; r</span>}}\end{array}\right.<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 8</span>}<span\; class="notranslate">\; )</span>$ Bring the linear combination of the extracted components Tk (k=1,...,r) from the HRV volume obtained in the first step analysis into the regression equation established by the M volume for r components, that is, t _r =ω _k1 HRV ₁ +...+ω _km HRV Substituting _m into M _j =t ₁ β _1j +...+t _r β _rj (j=1,...,p), the regression equation of standardized variables M _j =α _j1 HRV ₁ +...+α _jm HRV _m ; Finally, according to formula (3), L can be obtained. 3. the dual-source feature fusion ant colony tuning method of PID parameter in a kind of functional electric stimulation according to claim 1, it is characterized in that, described encoding is, according to knee joint angle and current pattern, and actual error situation Set the three parameters of PID to 5 significant figures, among which K _p has 2 digits before the decimal point and 3 digits after the decimal point; K _i and K _d have 1 digit before the decimal point and 4 digits after the decimal point. The specific coding is as follows: K _i = y _{6, j} × 10 ⁰ + y _{7, j} × 10 ^-1 + y _{8, j} × 10 ^-2 + y _{9, j} × 10 ^-3 + y _{10, j} × 10 ^-4 K _p = y _{1, j} × 10 ¹ + y _{2, j} × 10 ⁰ + y _{3, j} × 10 ^-1 + y _{4, j} × 10 ^-2 + y _{5, j} × 10 ^-3 K _d = y _{6, j} × 10 ⁰ + y _{7, j} × 10 ^-1 + y _{8, j} × 10 ^-2 + y _{9, j} × 10 ^-3 + y _{10, j} × 10 ^-4 And according to the output overshoot and other errors to design the adaptation function, the specific function is: fit＝α1×σ+β1×t+c×error (9) Among them, σ is the overshoot, t is the rise time, error is the relative error between the output joint angle and the preset joint angle, α1=0.1, β1=0.8, c=2. 4. the dual-source feature fusion ant colony tuning method of PID parameter in a kind of functional electric stimulation according to claim 1, it is characterized in that, the concrete workflow of ant colony algorithm tuning PID controller parameter is: Step1: parameter initialization, make time t=0 and number of cycles N _max =0, set the maximum number of cycles N _cmax , and place m ants at the starting point; Step2: Set the number of ants and the number of cycles; Step3: Ants search randomly. After a crawl is over, decide which feature variables are selected as the actual input variables, modify the taboo table pointer, that is, move the ant to a new element after selection, and move the element to the taboo of the ant individual table; Step4: Calculate the probability calculated by the ant individual according to the state transition probability formula, and select elements according to this probability; Step5: If the ant element has not been traversed, go to Step3, otherwise go to Step6; Step6: Update the pheromone concentration of the characteristic variable with high classification accuracy of pheromone concentration is enhanced, and will be selected with a greater probability in the next search; Step7: Finish the adjustment when it is satisfied, and the setting is over.

Images