CN101837164B - 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

Dual-source feature fusion ant colony tuning method for PID parameters in functional electrical stimulation

technical field

The invention relates to the field of equipment for limb rehabilitation by electrical pulse stimulation, in particular to a dual-source feature fusion ant colony tuning method for PID parameters in functional electrical stimulation.

Background technique

Functional Electrical Stimulation (FES) is a technology that stimulates limb motor muscles and peripheral nerves through current pulse sequences to effectively restore or reconstruct part of the motor function of paraplegic patients. According to statistics, due to the weak regeneration ability of the spinal cord, there is currently no effective treatment method that can directly repair the injury for paralyzed patients with spinal cord injury. The implementation of functional rehabilitation training is an effective measure. The number of paralyzed patients with spinal cord injury is increasing year by year, and functional rehabilitation training is an urgently needed technology. In the 1960s, Liberson successfully used electrical stimulation of the peroneal nerve for the first time to correct the gait of hemiplegic patients with foot drop, and created a new way of functional electrical stimulation for motor and sensory function rehabilitation. At present, FES has become an important rehabilitation treatment method for restoring or reconstructing part of the motor function of paraplegic patients. However, how to precisely control the trigger timing and pulse current intensity of FES to ensure that the effect of electrical stimulation can accurately complete the predetermined functional action is still the key to the technology of FES. According to statistics, there are still few researches on FES trigger control methods at present, and according to the deviation between the effect and the predetermined action, the closed-loop control is used to automatically adjust the FES stimulation intensity and timing parameters, thereby greatly improving the accuracy and stability of the FES system. Effective control methods are still being explored.

The handle reactions vector (HRV) is based on the process of standing and walking with the help of the walker. The utility provided by the walker to the patient can actually be divided into three clear and independent parts: forward and backward force propulsion, left and right In fact, this can also be understood as a new concept proposed by the patient for the additional mechanical appeal to the outside world in order to maintain his normal standing and walking, that is, the role of the patient on the walker during the process of standing and walking The combination is simplified as a concentrated load, which is represented by two mechanical vectors at the centroid of the cross-section at the midpoint of the handle, as shown in Figure 1. The resultant force of the vector components on the x, y, and z axes can respectively represent the patient’s ability to walk The obtained force propulsion, force balance and force support level of the machine. Wherein, the positive direction of the x-axis set by the defined coordinate system is the right direction of the patient, the positive direction of the y-axis is the forward direction of the patient, and the positive direction of the z-axis is the upward direction of the patient. In this way, the definition formula of HRV can also be written as:

[HRV] = [HRV _l , HRV _r ] ^T = [F _lx , F _ly , F _lz , F _rx , F _ry , F _rz ] ^T (1)

At present, HRV is widely used to monitor the patient's walking condition during electrical stimulation, and then prevent the patient from falling and causing secondary injury. This patent proposes to use this parameter to predict the knee joint angle, and then precisely control the current level of the FES system to ensure that the electrical stimulation effect can accurately complete the predetermined functional action and prevent muscle fatigue.

Proportional-integral-differential (PID) is a very practical feedback adjustment algorithm, which uses proportional, integral, and differential operations to obtain the required adjustment amount for feedback control of the system according to system detection or operating deviation. It is easy to operate and widely used in engineering practice. Especially when the characteristic parameters of the controlled system are not clear or it is difficult to measure online in time, the PID tuning algorithm can be used for safe closed-loop control. Facing the complex and time-varying operating environment of muscles, PID is still widely used in the field of functional electrical stimulation due to its good stability and reliable operation. The core technology of PID is to precisely determine the proportion, integral, and differential coefficients. Especially in the field of FES, the requirements for system stability are extremely strict, so the selection of PID parameters is particularly important. In order to achieve a better control effect in PID control, the three control functions of proportional, integral and differential must be adjusted to form a relationship of mutual cooperation and mutual restriction in the control quantity.

Contents of the invention

In order to overcome the shortcomings of the existing technology, a dual-source feature fusion ant colony tuning method for PID parameters in functional electrical stimulation is provided, which can accurately and stably control the current intensity of the FES system in real time, and effectively improve the accuracy and stability of the FES system , and obtain considerable social and economic benefits. In order to achieve the above object, the technical solution adopted in the present invention is: the dual-source feature fusion ant colony tuning method of PID parameters in the functional electrical stimulation, comprising 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.

The corresponding relationship between the actual joint angle and the 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, HRV0 and M0 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, and 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 X0 and Y0, repeat the extraction of components, and so on, assuming that r components are finally extracted, the regression equation of HRV0 and M0 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">\; \&Center\; Dot;</span><span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; \&CenterDot;</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">\; \&CenterDot;</span><span\; class="notranslate">\; \&Center\; Dot;</span><span\; class="notranslate">\; \&CenterDot;</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), that is, the regression equation of standardized variables M _j =α _j1 HRV ₁ +...+α _jm HRV _m ;

Finally, according to formula (3), L can be obtained.

The encoding is to set the three parameters of PID according to the knee joint angle and current mode, and the actual error situation as 5 significant figures, wherein K _p has 2 digits before the decimal point and 3 digits after the decimal point; K _i and K _d decimal points The first 1 digit, 4 digits after the decimal point, the specific encoding is as follows:

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 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;

The specific workflow of the ant colony algorithm to tune the parameters of the PID controller is as follows:

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.

The feature of the present invention is: use the HRV change of the walker to predict the change of the knee joint angle, then optimize the proportional coefficient, differential coefficient and integral coefficient of the PID through the ant colony algorithm, and then control the current pulse intensity of the FES system, effectively improving the FES system accuracy and stability.

Description of drawings

Figure 1 Schematic diagram of the definition of handle reaction vector (HRV).

Figure 2 is a block diagram of the HRV-based FES system.

Fig. 3 is a structural block diagram of ant colony algorithm tuning PID parameter control method.

Figure 4 Human model in functional electrical stimulation for walking aid.

Fig. 5 Schematic diagram of ant colony algorithm tuning PID parameter encoding.

Figure 6 Experimental scene.

Fig. 7 PID control tracking results of ant colony algorithm tuning.

Detailed ways

The gist of the present invention is to propose a new precision control method for FES, utilize the error of the knee joint angle predicted by the HRV parameter of the walker and the joint angle predicted by the actual knee joint angle, optimize the proportional coefficient of PID by ant colony algorithm, Integral coefficient and differential coefficient, and then accurately and stably control the current intensity of the FES system in real time. The invention can effectively improve the accuracy and stability of the FES system, and obtain considerable social and economic benefits.

The application structure of the new precision control technology in the HRV-based functional electrical stimulation walking aid is shown in Figure 2. Its workflow is as follows: firstly, use the HRV during the walking aid process to predict the knee joint angle, and secondly, use the ant colony algorithm to set PID parameters, real-time regulation of FES current level intensity. The schematic diagram of its tuning structure is shown in Figure 3, which is as follows: firstly, according to the upper and lower bounds of the value ranges of the three decision variables K _p , K _i and K _d of PID, determine the parameters such as the ant colony size and the dimension of the search space, and then It encodes, and then uses the corresponding relationship between the actual joint angle and the output joint angle of the muscle model as the fitness value calculated by the moderate evaluation function, and determines the parameter setting of the ant colony algorithm, and uses the ant random search to make its variables optimize the K of the PID The three coefficients _p , K _i and K _d use the fitness function to adjust the ant's search path each time and judge whether it reaches the preset goal. 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 the system output yout and its deviation from the muscle model under the new PID coefficient, and then enter the next step of ant colony algorithm self-learning and weight coefficient self-adjustment. Repeat this process, and finally realize the adaptive online tuning of PID control parameters, and use it in the FES system.

1 HRV prediction model of knee joint angle

During the walking aid process, when the user lifts his legs and takes a step under the action of functional electrical stimulation, in order to support the stability of the body, the force exerted by the user on the walker is different, because the different sizes of the joints will make the center of gravity of the human body In different positions, the force applied to overcome gravity is also different, and the plane position of the human body is also changed. In order to tilt the position, the force exerted on the plane also changes. Therefore, the joint angle and the user's perception of the walker The applied force has a certain relationship, as shown in Figure 4.

M＝L HRV+wPW (1)

where M is the knee joint angle, HRV is the handle reaction vector of the force exerted by the user on the walker, L is the relationship between HRV and M, w is the coefficient, W is the center of gravity of the upper arm, trunk and lower limbs, and P is the triple The relationship between heart and M.

In practice, due to the action of the walker, the movement of the center of gravity of the human body is small, and the angle of the knee joint can be expressed as

M＝L HRV (2)

Among them, M represents the knee joint angle, HRV represents the handle reaction vector of the force exerted by the user on the walker, and L represents the relationship between HRV and M. According to Equation 2, determining L can use HRV to obtain the knee joint angle at the corresponding moment.

L＝M×HRV ^-1 (3)

When the present invention solves L, has adopted the method for partial least square regression.

There are m HRV variables HRV1, . . . , HRVm, p M variables, M1 , . . . , Mp, and a data set of i (i=1, . T and U are components extracted from HRV variables and M variables respectively, and the components extracted here are usually called partial least squares factors.

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 model effect weight, and v ₁ =(v ₁₁ ,...,v _1p ) is the M variable weight. In order to ensure that T1 and U1 extract as much variation information as possible from the variable group they belong to, and at the same time ensure that the correlation between the two reaches the maximum, according to the property that the covariance of the components can be calculated by the inner product of the score vectors of the corresponding components, the above The requirement of extracting the first component is transformed into a conditional extremum problem.

得到。Among them, t1 and u1 are the score vectors of the first pair of components obtained from the samples, and HRV0 and M0 are the initial variables. 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, finding the matrix HRV′ ₀ M ₀ M′ ₀ HRV ₀ 's eigenvalue and eigenvector, its maximum eigenvalue is θ ₁ ² , the corresponding unit eigenvector is the solution ω ₁ , 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 first component extracted cannot reach the accuracy of the regression model, use the residual matrix E1 and F1 to replace X0 and Y0, repeat the extraction of components, and so on. Assuming that r components are finally extracted, the regression equation of HRV0 and M0 for r components is:

$\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">\; \&CenterDot;</span><span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; \&CenterDot;</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">\; \&Center\; Dot;</span><span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; \&CenterDot;</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 is obtained.

Finally, according to formula 3, L can be obtained.

The ant colony algorithm controls the PID parameters. First, encode the three parameters of the PID. According to the knee joint angle, current mode, and actual error, the three parameters of the PID are set to 5 significant figures, among which K _p before the decimal point 2 digits, 3 digits after the decimal point; K _i and K _d 1 digit before the decimal point, 4 digits after the decimal point. Its specific coding schematic diagram is shown in Fig. 5 .

as the 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 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.

The specific workflow of the ant colony algorithm to tune the parameters of the PID controller is as follows:

Step1: Parameter initialization. Set the time t=0 and the 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 in the 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 pheromone concentration The pheromone concentration of the characteristic variable with high classification accuracy is enhanced, and it will be selected with a greater probability in the next search

Step7: Finish the adjustment when it is satisfied, and the setting is over.

2. Control of PID parameter setting by ant colony algorithm

PID is composed of three parts: proportional unit P, integral unit I and differential unit D. According to the error of the system, the system is controlled by setting three parameters K _p , K _i and K _d .

$\mathrm{<span\; class="notranslate">\; yout</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}\mathrm{<span\; class="notranslate">\; error</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; t</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\mathrm{<span\; class="notranslate">\; error</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; error</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; error</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 9</span>}<span\; class="notranslate">\; )</span>$

Among them, K _p is the proportional coefficient, K _i is the integral coefficient, K _d is the differential coefficient, error is the deviation between the preset output and the actual output, u(t) is the output of PID, and it is also the input of the controlled system.

From the PID output formula (1), we can get

$\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}\mathrm{<span\; class="notranslate">\; error</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\mathrm{<span\; class="notranslate">\; error</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; error</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; error</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 10</span>}<span\; class="notranslate">\; )</span>$

according to:

Δu(t)=u(t)-u(t-1)

=K _p (error(t)-error(t-1))+K _i error(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 _i error(t)+K _d (error(t)-2error(t-1)+error(t-2 ))

……………………(12)

The invention adopts the ant colony algorithm to carry out self-adaptive optimization of the PID control parameters, takes the three parameters of the PID as a combination, and utilizes the ant colony algorithm to optimize to solve the combination problem. Ant colony algorithm is a new type of bionic algorithm derived from the biological world of nature. When using ant colony algorithm to solve optimization problems, the optimization problem is first transformed into the shortest path problem. Each ant starts from the initial node N ₀₀ , N ₀₁ ... N _0n , and walks sequentially through one of the sub-nodes of N ₁ , N ₂ ..., until the end point N _k0 , N _k1 ... N _kn Composition path (N _0t N _1t ... N _kt ), t ∈ [0, 1, 2 ... 9]. Its path can represent a binary feasible solution. Each time an ant visits a city, it has the following characteristics:

State transition rule: The state transition rule used by the ant colony algorithm is a random proportional rule based on the TSP problem, which gives the probability that ant k located in city i chooses to move to city j,

Where τ _ij (i, j) is the fitness of (i, j), η _ij (i, j) = (10-|y(i)-y(i) ^* |)/10, y(i) ant The value at i during the group search, y(i) ^* is the value at i during the last search. α is the relative importance of residual information, and β is the relative importance of expected value.

In the ant colony algorithm, the selection method is

Among them, q is a random number uniformly distributed on [0, 1], and q ₀ is a parameter variable on [0, 1].

Global update rules: The ant algorithm has different update algorithms. The global update principle adopted by the ant colony system only allows the ants with the global optimal solution to release pheromone. neighborhood of the best path.

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

Among them, ρ is the volatility coefficient of information number, and L _gb is the global optimal path found so far

Partial update information: each ant also updates the information trace in the process of establishing a solution

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

where γ ∈ [0, 1].

The ant colony algorithm controls the PID parameters. First, encode the three parameters of the PID. According to the knee joint angle, current mode, and actual error, the three parameters of the PID are set to 5 significant figures, among which K _p before the decimal point 2 digits, 3 digits after the decimal point; K _i and K _d 1 digit before the decimal point, 4 digits after the decimal point. Its specific coding schematic diagram is shown in Fig. 5 .

as the 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 specific workflow of the ant colony algorithm to tune the parameters of the PID controller is as follows:

Step1: Parameter initialization. Set the time t=0 and the 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 in the 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 pheromone concentration The pheromone concentration of the characteristic variable with high classification accuracy is enhanced, and it will be selected with a greater probability in the next search

Step7: Finish the adjustment when it is satisfied, and the setting is over.

3 experimental protocol

The experimental device uses a wireless transmission walker system and a Parastep functional electrical stimulation system produced by SIGMEDICS of the United States. The system includes a microprocessor and a stimulation pulse generation circuit, includes six stimulation channels, and is powered by batteries. The content of the experiment is: use the FES system to stimulate the relevant muscles of the lower limbs, so that the subjects follow the predetermined movements, and at the same time record the HRV applied on the walker. The voltage signal converted by the bridge network and the trajectory of the knee joint angle. The subjects are required to be in good health, without lower extremity muscle and bone diseases, without neurological diseases and severe cardiopulmonary diseases. During the experiment, the subjects sat in front of the walker,

Fix the stimulating electrodes to the corresponding parts, and keep the subjects relaxed when no electrical stimulation is applied. The FES experimental scene is shown in Fig. 5. The electrical stimulation pulse sequence adopts the classic Lilly waveform, the pulse frequency is 25Hz, the pulse width is 150μs, and the pulse current is adjustable within the range of 0-120m. In the experiment, HRV was recorded in real time to

And the stimulation intensity can be adjusted by changing the magnitude of the pulse current to change the angle of the knee joint generated by the stimulation. Before the experiment, set the expected trajectory of the knee joint angle, and use the goniometer to detect the change of the knee joint opening angle in real time during the experiment. The sampling rate of the experimental data is 128Hz, and the data recording time is 60s.

Beneficial effect

The new algorithm of ant colony algorithm to adjust PID parameters measures and adjusts the amplitude of FES pulse current, so that the knee joint angular motion generated by FES is close to the expected motion track. Figure 7 shows the tracking results of the PID control set by the ant colony algorithm. The red line in the figure represents the expected motion trajectory, and the blue line represents the actual output joint angle. The X-axis is time, and the Y-axis is the knee joint motion angle. In order to observe the control error of the ant colony algorithm tuning PID more clearly, as shown in Figure 8, the relative error between the preset input knee joint angle and the actual knee joint angle under the ant colony algorithm tuning PID control, it can be seen that the errors are all 5% Within, precise control can be achieved.

The gist of the present invention is to propose a new precision control method for FES, utilize the error of the knee joint angle predicted by the HRV parameter of the walker and the joint angle predicted by the actual knee joint angle, optimize the proportional coefficient of PID by ant colony algorithm, Integral coefficient and differential coefficient, and then accurately and stably control the current intensity of the FES system in real time. The invention can effectively improve the accuracy and stability of the FES system, and obtain considerable social and economic benefits. The best implementation plan is to use patent transfer, technical cooperation or product development.

Claims (3)

1. the double source feature fusion ant colony setting method of pid parameter in the functional electrostimulation is characterized in that, comprises the following steps:

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

Secondly, utilize ant group algorithm tuning PID parameter, real-time monitoring FES levels of current intensity, its flow process of adjusting is: at first according to three decision variable K of PID controller _p, K _iAnd K _dThe bound of span; Confirm to comprise the parameter of ant crowd population size, search volume dimension; And it is encoded; Utilize then through actual joint angles and muscle model HRV and export the fitness value of the corresponding relation of joint angles as the fitness evaluation function calculation, corresponding relation is L=M * HRV ^-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, adopts the method for PLS to confirm muscle model HRV output joint angles; Adopt ant group algorithm that the PID controller parameter is controlled, promptly confirm the parameter setting of ant group algorithm, utilize the ant random search to make its variable optimize three decision variable K of PID controller _p, K _iAnd K _d, utilize the fitness evaluation function adjusting each searching route of ant and judge whether to reach goal-selling, as reach goal-selling, calculate three decision variable K that final best position promptly gets the PID controller _p, K _iAnd K _d, as do not reach re-set target continuation optimizing, up to reaching goal-selling; Computing function property electro photoluminescence FES system output yout under the new PID controller parameter and with the deviation of muscle model HRV after get into the self study of next step ant group algorithm and three decision variable K of PID controller again _p, K _iAnd K _dThe weighting coefficient self-adjusting;

Last repeatedly step finally realizes the self-adaptation on-line tuning of PID controller parameter, and is used for functional electrostimulation FES system.

2. the double source feature fusion ant colony setting method of pid parameter is characterized in that in a kind of functional electrostimulation according to claim 1, adopts the method for PLS to find the solution L:

Be provided with m HRV variable HRV ₁..., HRV _m, p M variable, M ₁..., M _p, 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 concentrate from original variable and extract first couple of composition T ₁, U ₁Linear combination be:

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:

T wherein ₁, u ₁Be the score vector of first pair of composition of trying to achieve, HRV by sample ₀, M ₀Be 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 v1 is by formula

Obtain;

Next sets up the equation of initializaing variable to T1

$\left\{\begin{array}{c}{\mathrm{HRV}}_0=t_1{\mathrm{\&alpha;}}_1^{\&prime;}+E_1\\ M_0=t_1{\mathrm{\&beta;}}_1^{\&prime;}+F_1\end{array}\right.---\left(7\right)$

T wherein ₁Meaning is the same, α ' ₁=(α ₁₁..., α _1m), β ' ₁=(β ₁₁..., β _1p) be M variable t only ₁The time parameter vector, E ₁, F ₁Be 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;

Wherein, E ₁=HRV ₀- _t1 α ' ₁,

F ₁＝HRV ₀-t ₁β′ ₁

Can not reach the precision of regression model, utilization residual error battle array E like the first pair of composition that extracts ₁, F ₁Replace HRV ₀, M ₀, repeat to extract composition, and the like, suppose finally to have extracted r composition, HRV ₀, M ₀Regression equation to r composition is:

$\left\{\begin{array}{c}{\mathrm{HRV}}_0=t_1{\mathrm{\&alpha;}}_1^{\&prime;}+...+t_r{\mathrm{\&alpha;}}_r^{\&prime;}+E_r\\ M_0=t_1{\mathrm{\&beta;}}_1^{\&prime;}+...+t_r{\mathrm{\&beta;}}_1^{\&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 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

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

3. the double source feature fusion ant colony setting method of pid parameter is characterized in that in a kind of functional electrostimulation according to claim 1, and 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 ant placed starting point;

Step2: ant number and cycle index are set;

Step3: ant 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 ant is moved to new element, and move to this element in the individual taboo table of ant;

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

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

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

Step7: satisfy termination condition, the end of adjusting.

Images