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

Abstract

The invention relates to the field of equipment for limb rehabilitation by electric pulse stimulation. Provide a dual-source feature fusion chaotic particle swarm tuning method for PID parameters in functional electrical stimulation, 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. The technical solution adopted by the invention Yes, firstly, use the handle reaction vector HRV of the walking aid process to predict the knee joint angle; secondly, use the chaotic particle swarm optimization algorithm to set the proportional calculus PID parameters, adjust the FES current level and intensity in real time, and finally realize the adaptive proportional calculus PID control parameters Online tuning, and used for functional electrical stimulation FES system. The invention is mainly applied to setting PID parameters in functional electrical stimulation.

Description

Functional Electrical Stimulation PID Parameters Dual-source Feature Fusion Particle Swarm Tuning Method

technical field

The invention relates to the field of equipment for limb rehabilitation by means of electrical pulse stimulation, in particular to a method for setting dual-source feature fusion chaotic particle groups of 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 is still little research on the trigger control of FES 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 ₁ , 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 deficiencies of the existing technology, a dual-source feature fusion chaotic particle swarm 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 chaotic microparticle swarm tuning method of PID parameters in functional electrical stimulation, including:

First, the knee joint angle is predicted using the handle reaction vector HRV of the walking aid process;

Secondly, use the chaotic particle swarm optimization algorithm to tune the proportional calculus PID parameters, and adjust the FES current level in real time. The tuning process is as follows: first, according to the upper and lower bounds of the value ranges of the three decision variables K _p , K _i and K _d of the proportional calculus PID , determine the parameters such as the size of the particle swarm, the dimension of the search space, and initialize the velocity and position of the particle swarm, and then use the corresponding relationship between the actual joint angle and the output joint angle of the muscle model as a moderate evaluation function to calculate the value of each particle in the particle swarm The fitness value, and compare its fitness with its own best position fitness value, and use it as the representative value of the particle, and then adjust the speed and other parameters of the particle to change the best position of the particle until it is stable. Calculate the final optimal position to get the three coefficients K _p , K _i and K _d of the proportional calculus PID, calculate the system output yout and its deviation from the muscle model output joint angle under the new proportional calculus PID coefficients Enter the next step of neural network self-learning and weighting coefficient self-adjustment, repeat this process, and finally realize the adaptive online tuning of proportional calculus PID control parameters, and use it in the functional electrical stimulation FES system.

The muscle model output joint angle adopts the method of partial least squares regression, namely:

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 variables, 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 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 ₀ , the largest eigenvalue of which is θ ₁ ² , and the corresponding unit eigenvector is the obtained 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">\; <figure-callout\; id="0"\; label="HRV"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">HRV</figure-callout></span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; <figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout></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">\; <figure-callout\; id="0"\; label="m"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">m</figure-callout></span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; <figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout></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 t ₁ has the same meaning as before, α ₁ ′=(α ₁₁ ,…,α _1m ), β ₁ ′=(β ₁₁ ,…,β _1p ) is the parameter vector when there is only one M quantity t1, E1 and F1 are respectively n×m and n×p residual matrix, according to the ordinary least squares method, the coefficient vectors α ₁ and β ₁ can be obtained, 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">\; <figure-callout\; id="0"\; label="HRV"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">HRV</figure-callout></span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; <figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout></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">\; <figure-callout\; id="0"\; label="E.\; r\; m"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">E.</figure-callout></span><figure-callout\; id="0"\; label="E.\; r\; m"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">\; </figure-callout>}}_{\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">\; <figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout></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 the formula L=M×HRV ^-1 , L can be obtained. 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.

The PID parameter tuning by using the chaotic particle swarm optimization algorithm is further refined as:

Proportional calculus 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">\; [</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">\; ]</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. output formula $\mathrm{yout}\left(t\right)=K_p\mathrm{error}\left(t\right)+K_i\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}\mathrm{error}\left(j\right)+K_d[\mathrm{error}\left(t\right)-\mathrm{error}(t-1)]$ can 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">\; [</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">\; ]</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 chaotic particle swarm optimization algorithm is used for the adaptive optimization of proportional calculus PID control parameters. The search space is selected as 3-dimensional, that is, the three parameters of the PID controller. The group size m=20 is selected. Randomly generated within the space range of , respectively expressed as: v _i = (v _i1 , v _i2 , v _i3 ), x _i = (x _i1 , x _i2 , x _i3 ), remember that the i-th particle has searched for the optimal position so far is p _i =(p _i1 , p _i2 , p _i3 ), the optimal position searched so far by the entire particle swarm is p _gi =(p _gi1 , p _gi2 , p _gi3 ), where i=1, 2,... , 20, the particle swarm optimization algorithm uses the following formula to operate on the particle swarm,

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

x _id ← x _id +v _id (14)

Among them, i=1, 2, ..., 20; the learning factors c ₁ , c ₂ and c ₃ are non-negative numbers, generally set at 0.5; r ₁ , r ₂ and r ₃ are between [0, 1] Random number, q _id is the position of randomly selected particles;

The specific steps to realize are:

1. Determine the parameters: learning factors c ₁ , c ₂ and c ₃ , and the size N of the group, the times of evolution and the times of chaos optimization;

2. Randomly generate N particles for operation;

3. Operate the particles according to formulas (13) and (14);

4. Perform chaos optimization on the optimal position p _gi =(p _gi1 , p _gi2 , p _gi3 ), and map p _gid (i=1, 2, ..., 20) to

m＝1，2，…再把产生的混沌变量序列

(m＝1，2，…)通过逆映射

(m＝1，2，…)返回到原解空间，得

(m＝1，2，…)The definition domain [0, 1] of Logistic equation z _i+1 = μ z _i (1-z _i ) i = 0, 1, 2, ...; i=1, 2,..., 20, then, use the Logistic equation to iterate to generate chaotic variables

m=1, 2, ... and then the generated chaotic variable sequence

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

(m=1, 2, ...) Return to the original solution space, get

(m=1, 2,...)

(m＝1，2，…)计算其适应值，保留性能最好的可行解p^*；Every feasible solution experienced by the chaotic variable in the original solution space

(m=1, 2, ...) calculates its fitness value, retains the best feasible solution p ^* of performance;

5. A particle randomly selected from the current population is replaced by p ^* ;

6. If the maximum algebra is reached or a satisfactory solution is obtained, the optimization process ends, otherwise, return to step 3.

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 chaotic particle swarm 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 the control method of PID parameter tuning by chaotic particle swarm algorithm.

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

Figure 5 Experimental scene.

Figure 6 PID control tracking results of chaotic particle swarm optimization algorithm tuning.

Figure 7 The relative error between the preset input closing angle and the actual output under the control of the PID parameters adjusted by the chaotic particle swarm algorithm.

Detailed ways

The application structure of the new precision control technology in the HRV-based functional electrical stimulation walking aid is shown in Figure 2, and its workflow is as follows: firstly, the knee joint angle is predicted using the HRV during the walking aiding process, and secondly, the chaotic particle swarm optimization algorithm is used Adjust the PID parameters and adjust the FES current level in real time. 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 particle swarm group size and the dimension of the search space, and initialize The velocity and position of the particle swarm, and then use the corresponding relationship between the actual joint angle and the output joint angle of the muscle model as a moderate evaluation function to calculate the fitness value of each particle in the particle swarm, and adapt its fitness to its own best position Then, adjust the particle speed and other parameters, change the best position of the particle until it is stable, and calculate the final best position to get K _p , K _i and PID K _d three coefficients. Calculate the system output yout and its deviation from the muscle model under the new PID coefficient, and then enter the next step of neural network 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 raises 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 solving L in this patent, the method of partial least squares regression is adopted.

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 the unit vectors ω ₁ and v ₁ , so that θ ₁ = ω ₁ HRV ₀ ′M ₀ v ₁ is the largest, that is, finding the matrix HRV ₀ ′M ₀ M ₀ ′HRV ₀ The eigenvalues and eigenvectors of , 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">\; <figure-callout\; id="0"\; label="HRV"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">HRV</figure-callout></span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; <figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout></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">\; <figure-callout\; id="0"\; label="m"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">m</figure-callout></span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; <figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout></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 respectively n ×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">\; <figure-callout\; id="0"\; label="HRV"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">HRV</figure-callout></span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; <figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout></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">\; <figure-callout\; id="0"\; label="E.\; r\; m"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">E.</figure-callout></span><figure-callout\; id="0"\; label="E.\; r\; m"\; filenames="HDA0000021781120000031.png"\; state="\{\{state\}\}">\; </figure-callout>}}_{\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">\; <figure-callout\; id="1"\; label="t"\; filenames="HDA0000021781120000032.png"\; state="\{\{state\}\}">t</figure-callout></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 is obtained.

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

2. Control of PID parameter tuning by chaotic particle swarm 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">\; [</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">\; ]</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">\; [</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">\; ]</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 present invention adopts chaotic particle swarm algorithm to carry out self-adaptive optimization of PID control parameters, selects the search space to be 3 dimensions, namely respectively the three parameters of the PID controller, selects the group scale m=20, and the initial speed and position of the group are in a certain Randomly generated within the space range, respectively expressed as: v _i = (v _i1 , v _i2 , v _i3 ), x _i = (x _i1 , x _i2 , x _i3 ), record the optimal position of the i-th particle searched so far as p _i =(p _i1 , p _i2 , p _i3 ), the optimal position searched so far by the entire particle swarm is p _gi =(p _gi 1 , p _gi2 , p _gi3 ). Among them, i=1, 2, ..., 20, the particle swarm optimization algorithm uses the following formula to operate on the particle swarm.

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

x _id ← x _id +v _id (14)

Among them, i=1, 2, ..., 20; the learning factors c ₁ , c ₂ and c ₃ are non-negative numbers, generally set at 0.5; r ₁ , r ₂ and r ₃ are between [0, 1] Random number, q _id is the position of randomly selected particles.

The specific steps for its realization are:

1. Determine the parameters: learning factors c ₁ , c ₂ and c ₃ , and the size N of the group, the times of evolution and the times of chaos optimization;

2. Randomly generate N particles for operation;

3. Operate the particles according to formulas (13) and (14);

i＝1，2，…，20，然后，用Logistic方程进行迭代产生混沌变量m＝1，2，…，再把产生的混沌变量序列(m＝1，2，…)通过逆映射

(m＝1，2，…)返回到原解空间，得

(m＝1，2，…)4. Perform chaos optimization on the optimal position p _gi = (p _gi1 , p _gi2 , p _gi3 ), and map p _gid (i=1, 2,..., 20) to the Logistic equation z _i+1 = μ z _i ( 1-z _i ) the domain [0, 1] of i=0, 1, 2, ...;

i=1, 2,..., 20, then, use the Logistic equation to iterate to generate chaotic variables m=1, 2, ..., and then the generated chaotic variable sequence (m=1, 2, ...) by inverse mapping

(m=1, 2, ...) Return to the original solution space, get

(m=1, 2,...)

(m＝1，2，…)计算其适应值，保留性能最好的可行解p^*。Every feasible solution experienced by the chaotic variable in the original solution space

(m=1, 2, . . . ) calculate its fitness value, and keep the feasible solution p ^* with the best performance.

5. A particle randomly selected from the current population is replaced by p ^* .

6. If the maximum algebra is reached or a satisfactory solution is obtained, the optimization process ends, otherwise return to step 3.

3. Experimental plan

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, and the stimulating electrodes were fixed on the corresponding parts. When no electrical stimulation was applied, the subjects remained relaxed. 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 is recorded in real time and the stimulation intensity can be adjusted by changing the magnitude of the pulse current to change the knee joint angle 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 chaotic particle swarm algorithm to tune 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 trajectory. Figure 6 is the tracking result of PID control tuned by chaotic particle swarm optimization 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 PID tuning by chaotic particle swarm optimization algorithm more clearly, as shown in Figure 7, the relative error between the preset input knee joint angle and the actual knee joint angle under the control of chaotic particle swarm optimization PID, it can be seen that the errors are all within 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, and optimize the proportional coefficient of PID through the chaotic particle swarm 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 (2)

1. A functional electric stimulation PID parameter dual-source feature fusion microparticle swarm tuning method is characterized in that it comprises the following steps: First, the knee joint angle is predicted using the handle reaction vector HRV of the walking aid process; _Secondly , use the chaotic particle swarm optimization algorithm to adjust the _proportional calculus PID parameters, and adjust the _FES current level in real time. Boundary, determine the parameters including the size of the particle swarm and the dimension of the search space, and initialize the velocity and position of the particle swarm, and then use the corresponding relationship between the actual joint angle and the output joint angle of the muscle model as a moderate evaluation function to calculate each of the particle swarm The fitness value of the particle, the corresponding relationship is L=M×HRV ^-1 , where 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, The partial least squares regression method is used to determine the output joint angle of the muscle model, and the fitness value of the particle is compared with the best position fitness value, and the fitness value of the particle and the best position of the particle itself are taken as The representative value of the particle, and then adjust the speed of the particle and other parameters, change the optimal position of the particle until it is stable, and calculate the final optimal position to get the three decision variables K _p , K _i and K of the proportional calculus PID parameters _d , under the new proportional calculus PID parameters, calculate the output yout of the functional electrical stimulation FES system and its deviation from the output joint angle of the muscle model, and then enter the next step of the self-learning of the chaotic particle swarm algorithm neural network and the three-dimensional PID parameters at that time The self-adjustment of the weighting coefficients of the decision variables K _p , K _i and K _d is repeated, and finally the adaptive online tuning of the proportional calculus PID parameters is realized and used in the functional electrical stimulation FES system. 2. a kind of functional electrical stimulation PID parameter dual-source feature fusion microparticle swarm tuning method according to claim 1, is characterized in that, muscle model output joint angle is the method that adopts partial least squares regression, namely: There are m HRV variables HRV1, ..., HRVm, p M variables, M1, ..., Mp, a data set of i observations in total, i=1, ..., n; T and U are respectively from HRV variables and M The components extracted from the variables are called partial least squares factors, and 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, and the above-mentioned requirement for extracting the first pair of components is transformed into the condition Extreme value problem:

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 ₀ , the largest eigenvalue of which is θ ₁ ² , and the corresponding unit eigenvector is The desired solution ω ₁ , and v ₁ is given by the formula get; Secondly, establish the equation of the initial variable to T1 Where t1 has the same meaning as before, α′ ₁ =(α ₁₁ ,…,α _1m ), β′ ₁ =(β ₁₁ ,…,β _1p ) is the parameter vector when there is only one M variable t1, and E1 and F1 are n respectively ×m and n×p residual matrix, the coefficient vectors α ₁ and β ₁ are obtained according to the ordinary least square method, where α ₁ is called the model effect load; If the first pair of components extracted cannot reach the accuracy of the regression model, use the residual matrix E1, F1 instead of HRV ₀ , M ₀ , repeat the extraction of components, and so on, assuming that r components are finally extracted, HRV0, M0 pair r components The regression equation for is:

Bring the linear combination of the extracted components Tk (k=1,...,r) into the regression equation established by M variables for r components from the HRV volume obtained in the first step analysis, 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, L can be obtained according to the formula L=M×HRV ^-1 . 3. A kind of functional electrical stimulation PID parameter dual-source feature fusion microparticle swarm tuning method according to claim 1, is characterized in that, utilizes chaotic particle swarm algorithm tuning PID parameter, further refines as: Proportional calculus PID is composed of three parts: proportional unit P, integral unit I and differential unit D. According to the error of the functional electrical stimulation system, the functional electrical stimulation system can be adjusted by setting three parameters K _p , K _i and K _d Take control:

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. output formula

can get

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 chaotic particle swarm optimization algorithm is used for the adaptive optimization of the proportional calculus PID control parameters. The search space is selected as 3-dimensional, that is, the three parameters of the PID controller, and the group size m=20 is selected. Randomly generated within the space range, respectively expressed as: v _i = (v _i1 , v _i2 , v _i3 ), x _i = (x _i1 , x _i2 , x _i3 ), record the optimal position of the i-th particle searched so far as p _i =(p _i1 , p _i2 , p _i3 ), the optimal position searched so far by the entire particle swarm is p _gi =(p _gi1 , p _gi2 , p _gi3 ), where i=1, 2,..., 20. The particle swarm optimization algorithm uses the following formula to operate particle swarms, v _id ←v _id +c ₁ r ₁ (p _id -x _id )+c ₂ r ₂ (p _gid -x _id )+c ₃ r ₃ (q _id -x _id ) (13) x _id ← x _id +v _id (14) Among them, the learning factors c ₁ , c ₂ and c ₃ are non-negative numbers with a value of 0.5; r ₁ , r ₂ and r ₃ are random numbers between [0, 1], and q _id is the position of randomly selected particles , d represents the encoding method of ions, d=10; The specific steps to realize are: 1. Determine the parameters: learning factors c ₁ , c ₂ and c ₃ , and the size N of the group, the times of evolution and the times of chaos optimization; 2. Randomly generate N particles for operation; 3. Operate the particles according to formulas (13) and (14); 4. Perform chaos optimization on the optimal position p _gi = (p _gi1 , p _gi2 , p _gi3 ), map p _gid to the Logistic equation z _i+1 = μ _{z i} (1-z _i ), z _i+1 is chaos The value of the variable, its definition domain is [0, 1], μ is the Logistic control parameter, which controls the chaotic state;

Then, use the Logistic equation to iterate to generate chaotic variables

m=1, 2, ..., and then the generated chaotic variable sequence

by inverse mapping

Returning to the original solution space, we get

a _id represents each smallest parameter among PID parameters, b _id represents each largest parameter among PID parameters; Every feasible solution experienced by the chaotic variable in the original solution space

Calculate its fitness value and keep the feasible solution p ^* with the best performance; 5. A particle randomly selected from the current population is replaced by p ^* ; 6. If the maximum number of iterations is reached or a satisfactory solution is obtained, the optimization process ends, otherwise, return to step 3. the

Images