CN102629847A - Asynchronous motor pure electronic speed feedback method - Google Patents

Abstract

The invention provides an asynchronous motor pure electronic speed feedback method, which specifically comprises the following operation steps of: firstly obtaining a discrete time form expression of an asynchronous motor mathematical model, obtaining an asynchronous motor reduced order EKF (Extended Kalman Filter) speed estimation algorithm according to an extended Kalman filter algorithm, designing a FPGA (Field Programmable Gate Array) to realize a reduced order EKF speed estimation algorithm structure and carrying out hardware language VHDL (Verilog Hardware Description Language) description on the FPGA based on the described algorithm structure to obtain, send and transmit back state estimated values i alpha s, i beta s, Psi alpha r, Psi beta r and omega r to a main control DSP (Digital Signal Processor) through a port. The asynchronous motor pure electronic speed feedback method has the beneficial effects that the pressure of the main control DSP on real-time operation amount is greatly reduced, so that more storage spaces and more operation spaces are left for speed and current control; the EKF speed estimation algorithm realized by the FPGA in parallel is used for finishing within 1 microsecond, so that a less sampling period can be selected by using the speed estimation algorithm, therefore, the speed estimation precision is greatly improved.

Description

The pure electronic type speed feedback of asynchronous machine method

Technical field

The invention belongs to industrial drive technology field, relate to a kind of speed feedback method of pure electronic type, be specifically related to a kind of method that realizes the online estimation of asynchronous machine rotating speed with pure hardware circuit.

Background technology

In asynchronous machine high-performance speed regulating control, need feedback speed signal to do speed closed loop control, to obtain Fast Dynamic response and higher steady precision.

In recent years, power electronic technology and computer technology develop rapidly to realizing that the high-performance AC governing system provides necessary condition.Since rotating speed, magnetic flux transducer be installed increased the price of drive system, reduced system reliability, destroyed intrinsic robustness and the simplicity of asynchronous motor.Therefore, Speedless sensor control not only becomes a research direction of modern AC drive control technology, also is simultaneously the key technology of development high performance universal frequency-variable controller.The key problem of Speedless sensor control technology is that motor speed is accurately estimated, and gives speed control with the speed feedback of estimating.

Along with the develop rapidly of high-performance digital signal processor, various method for estimating rotating speed emerge in an endless stream.Generally speaking, can be divided into following several types: direct computing method; Directly synthetic by state equation; MRAS (ModelReference Adaptive System-model reference adaptive); Method based on the full scalariform attitude of self adaptation observer; Method based on EKF (Extended Kalman Filter-extended Kalman filter); Based on neural network method; The teeth groove Harmonic detection; High-frequency signal injection method etc.

Direct computing method obtains the expression formula of slip or rotating speed based on asynchronous machine Park equation from the electromagnetic relationship formula of motor.The direct computing method that just is based on the rotor back electromotive force that adopts in the Speedless sensor vector control product of early stage Toshiba, Hitachi, Fuji.

Adopt the direct computing method Calculation Speed, directly perceived simple, real-time is good.But, because these class methods belong to the open loop observation procedure in essence, there is not the error correction link, shortcoming is also fairly obvious.The problem that exists has: 1) owing to can not guarantee whether vector control can correctly realize in the dynamic process, thereby can not guarantee the dynamic property of whole drive system; 2) interference rejection ability of system is very poor, and steady-state behaviour receives the influence of load torque, parameter of electric machine disturbance bigger; 3) to method based on the rotor back electromotive force because the value of back-emf is very little during low speed, cause calculating speed error bigger.

Utilize the state equation of asynchronous machine to set up flux observer, can observe the rotor flux or the stator magnetic linkage of motor, just can obtain synchronous speed signal ω through the mode of differentiate _eIntroduce slip ω again _SlCalculating, through ω _r=ω _e-ω _SlObtain rotor electrical angle speed omega _rABB, Japan's motor once used this method in product.Obviously, directly synthetic method also can be very responsive to the parameter of electric machine, and owing to belong to the open loop observation procedure equally, so estimated performance is relatively poor.

The basic thought of estimating based on the rotating speed of model reference adaptive method is with the equation that does not contain rotating speed model as a reference; And the equation that will contain rotating speed is as adjustable model; Two models have the output of same physical meaning, utilize the suitable adaptive law of error structure of these two model output variables, the rotating speed in the real-time regulated adjustable model; With the output of the output tracking reference model that reaches adjustable model, finally realize the purpose that rotating speed is estimated.Model reference adaptive method is a kind of method based on stability Design, can guarantee the asymptotic convergence of estimating.According to the different choice of reference model and adjustable model, multiple MRAS method for estimating rotating speed is arranged.Use maximum MRAS rotating speed algorithm for estimating and be with the voltage model that is used for rotor flux observation model as a reference, with the current model that contains rotary speed information algorithm as adjustable model.Can make full use of the result of calculation of flux observation process like this, the realization that rotating speed is estimated only need increase very little amount of calculation.Because the working voltage model calculates rotor flux, has introduced pure integral element, makes the magnetic linkage model receive the influence of initial value for integral and drift serious, the rotating speed estimated result is inaccurate.And poor-performing during low speed.

Above-mentioned several kinds of method for estimating rotating speed all need use open loop reconstruct (observation) value of rotor flux, so all belong to the observation procedure of open loop character in essence.These open loop method for estimating rotating speed are to overcome the parameter of electric machine changes and integrator drift brings error, so estimated performance is relatively poor.Can't support its commercialization based on present level of hardware such as neural network method, teeth groove Harmonic detection, high-frequency signal injection methods.Therefore, repeat no more here.

At present, the rotating speed estimated performance should be the observation procedure of closed loop class preferably, like full order observer, extended Kalman filter etc.Full order observer utilizes the measured value of rotor flux measured value and stator current and deviation of measuring value to estimate rotating speed, according to the dynamical equation and the Lyapunov Theory of Stability derivation adaptive law of state error.The full scalariform attitude of common self adaptation observer has based on ELO's (Extended LuenbergerObserver-expansion Long Beige observer) and based on sliding mode observer.These class methods still contain the thought of MRAS, and just the reference model has here become asynchronous machine itself, and adjustable model has become the closed loop full order observer.These class methods are better in the stability and the dynamic characteristic in full speed degree territory.But ELO is based on the observation procedure of certainty equation, the influence that the rotating speed estimation effect is changed by the parameters of electric machine such as electric machine non-linear characteristic and rotor resistance still.Self adaptation method for estimating rotating speed based on sliding mode observer has good stability to parameter error, system noise, has improved the robustness of estimating.

Since nineteen sixty R.E.Kalman proposed KF (Kalman Filter-Kalman filter), it had obtained using widely in every field, has become the canonical algorithm of linear system state estimation problem.Aspect non-linear estimations; Carry out the EKF (extended Kalman filter) of Filtering Estimation through the first-order linear method with Kalman filter, the non-linear system status algorithm for estimating that rely on that method is simple, realizability strong, advantages such as convergence becomes extensive use the most fast.EKF provides a kind of non-linear estimations algorithm of iteration form, has avoided differentiating, and can come the convergence rate of adjustment state estimation through the regulating error covariance matrix.In addition, of paramount importancely be, owing to be to be based upon on the random process model of system based on the method for estimating rotating speed of EKF, therefore have very strong anti-noise ability, this is that other rotating speed algorithm for estimating is not available.The statistics essence of its algorithm makes EKF be very suitable for overcoming the uncertainty and the non-linearity of asynchronous machine model, and estimated performance is superior.The deficiency of EKF method for estimating rotating speed is that amount of calculation is unfavorable for canbe used on line more greatly.

Summary of the invention

The purpose of this invention is to provide the pure electronic type speed feedback of a kind of asynchronous machine method; Neither take master control DSP the memory space resource and running time resource, can in full speed range, provide in real time again than the model reference adaptive rotating speed and estimate more accurate speed feedback value.

The technical scheme that the present invention adopted is, the pure electronic type speed feedback of asynchronous machine method, and the concrete operations step is following:

The first step obtains the discrete time formal representation of asynchronous machine Mathematical Modeling:

At first according to two mutually static α β shaft models of asynchronous machine:

$\left[\begin{array}{c}{\stackrel{\·}{i}}_{\mathrm{\αs}}\\ {\stackrel{\·}{i}}_{\mathrm{\βs}}\\ {\stackrel{\·}{\mathrm{\Ψ}}}_{\mathrm{\αr}}\\ {\stackrel{\·}{\mathrm{\Ψ}}}_{\mathrm{\βr}}\end{array}\right]=\left[\begin{array}{cccc}-(\frac{R_s}{\mathrm{\σ}{L}_s}+\frac{1-\mathrm{\σ}}{\mathrm{\σ}{\mathrm{\τ}}_r})& 0& \frac{L_m}{\mathrm{\σ}{L}_s{L}_r{\mathrm{\τ}}_r}& {\mathrm{\ω}}_r\frac{L_m}{\mathrm{\σ}{L}_s{L}_r}\\ 0& -(\frac{R_s}{\mathrm{\σ}{L}_s}+\frac{1-\mathrm{\σ}}{\mathrm{\σ}{\mathrm{\τ}}_r})& -{\mathrm{\ω}}_r\frac{L_m}{\mathrm{\σ}{L}_s{L}_r}& \frac{L_m}{\mathrm{\σ}{L}_s{L}_r{\mathrm{\τ}}_r}\\ \frac{L_m}{{\mathrm{\τ}}_r}& 0& -\frac{1}{{\mathrm{\τ}}_r}& -{\mathrm{\ω}}_r\\ 0& \frac{L_m}{{\mathrm{\τ}}_r}& {\mathrm{\ω}}_r& -\frac{1}{{\mathrm{\τ}}_r}\end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{\αs}}\\ i_{\mathrm{\βs}}\\ {\mathrm{\Ψ}}_{\mathrm{\αr}}\\ {\mathrm{\Ψ}}_{\mathrm{\βr}}\end{array}\right]+$

$\left[\begin{array}{cc}\frac{1}{\mathrm{\σ}{L}_s}& 0\\ 0& \frac{1}{\mathrm{\σ}{L}_s}\\ 0& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{u}_{\mathrm{\αs}}\\ u_{\mathrm{\βs}}\end{array}\right]---\left(1\right)$

In the formula (1), R _s, R _rBe stator and rotor resistance parameters; L _s, L _rBe the stator and rotor inductance; L _mBe mutual inductance; σ is a magnetic leakage factor, σ=1-(L _m ²/ L _sL _r); τ _rBe rotor time constant, τ _r=L _r/ R _rω _rBe rotor electrical angle speed; u _{α s}, u _{β s}Be α, β axle stator voltage; i _{α s}, i _{β s}Be α, β axle stator current; Ψ _{α r}, Ψ _{β r}Be α, β axle rotor flux;

Differential equation group form with formula (1) writing formula (2)～(5):

${\stackrel{\·}{i}}_{\mathrm{\αs}}=-(\frac{R_s}{\mathrm{\σ}{L}_s}+\frac{1-\mathrm{\σ}}{\mathrm{\σ}{\mathrm{\τ}}_r})i_{\mathrm{\αs}}+\frac{L_m}{\mathrm{\σ}{L}_s{L}_r{\mathrm{\τ}}_r}{\mathrm{\Ψ}}_{\mathrm{\αr}}+{\mathrm{\ω}}_r\frac{L_m}{\mathrm{\σ}{L}_s{L}_r}{\mathrm{\Ψ}}_{\mathrm{\βr}}+\frac{1}{\mathrm{\σ}{L}_s}{u}_{\mathrm{\αs}}---\left(2\right)$

${\stackrel{\·}{i}}_{\mathrm{\βs}}=-(\frac{R_s}{\mathrm{\σ}{L}_s}+\frac{1-\mathrm{\σ}}{\mathrm{\σ}{\mathrm{\τ}}_r})i_{\mathrm{\βs}}-{\mathrm{\ω}}_r\frac{L_m}{\mathrm{\σ}{L}_s{L}_r}{\mathrm{\Ψ}}_{\mathrm{\αr}}+\frac{L_m}{\mathrm{\σ}{L}_s{L}_r{\mathrm{\τ}}_r}{\mathrm{\Ψ}}_{\mathrm{\βr}}+\frac{1}{\mathrm{\σ}{L}_s}{u}_{\mathrm{\βs}}---\left(3\right)$

${\stackrel{\·}{\mathrm{\Ψ}}}_{\mathrm{\αr}}=\frac{L_m}{{\mathrm{\τ}}_r}{i}_{\mathrm{\αs}}-\frac{1}{{\mathrm{\τ}}_r}{\mathrm{\Ψ}}_{\mathrm{\αr}}-{\mathrm{\ω}}_r{\mathrm{\Ψ}}_{\mathrm{\βr}}---\left(4\right)$

${\stackrel{\·}{\mathrm{\Ψ}}}_{\mathrm{\βr}}=\frac{L_m}{{\mathrm{\τ}}_r}{i}_{\mathrm{\βs}}+{\mathrm{\ω}}_r{\mathrm{\Ψ}}_{\mathrm{\αr}}-\frac{1}{{\mathrm{\τ}}_r}{\mathrm{\Ψ}}_{\mathrm{\βr}}---\left(5\right)$

Formula (2), (3) both sides are with taking advantage of σ L _s,

$\mathrm{\σ}{L}_s{\stackrel{\·}{i}}_{\mathrm{\αs}}=-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\τ}}_r})i_{\mathrm{\αs}}+\frac{L_m}{L_r{\mathrm{\τ}}_r}{\mathrm{\Ψ}}_{\mathrm{\αr}}+{\mathrm{\ω}}_r\frac{L_m}{L_r}{\mathrm{\Ψ}}_{\mathrm{\βr}}+u_{\mathrm{\αs}}---\left(6\right)$

$\mathrm{\σ}{L}_s{\stackrel{\·}{i}}_{\mathrm{\βs}}=-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\τ}}_r})i_{\mathrm{\βs}}-{\mathrm{\ω}}_r\frac{L_m}{L_r}{\mathrm{\Ψ}}_{\mathrm{\αr}}+\frac{L_m}{L_r{\mathrm{\τ}}_r}{\mathrm{\Ψ}}_{\mathrm{\βr}}+u_{\mathrm{\βs}}---\left(7\right)$

Suppose that change in rotational speed is very little in a sampling period, can use approx

Replace mechanical equation, reelect x=[Ψ _{α r}, Ψ _{β r}, ω _r] ^TMake state, u=[i _{α s}, i _{β s}] ^TControl, the transposition computing of subscript T representing matrix is deformed into output equation with formula (6), formula (7) on this basis; After the arrangement, the 3 scalariform attitude equations that obtain the asynchronous machine system are:

$\stackrel{\·}{x}=f(x,u)---\left(8\right)$

Promptly

$\left[\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&alpha;r}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&beta;r}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_r\end{array}\right]=\left[\begin{array}{ccc}-\frac{1}{{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}& 0\\ {\mathrm{\&omega;}}_r& -\frac{1}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}\\ {\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}\\ {\mathrm{\&omega;}}_r\end{array}\right]+\left[\begin{array}{cc}\frac{L_m}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\\ 0& \frac{L_m}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}\\ 0& 0\end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{\&alpha;s}}\\ i_{\mathrm{\&beta;s}}\end{array}\right]---\left(9\right)$

The measurement equation is:

y＝h(x) (10)

Promptly

$y=\left[\begin{array}{c}{u}_{\mathrm{\&alpha;s}}-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}\\ u_{\mathrm{\&beta;s}}-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}\end{array}\right]=\left[\begin{array}{ccc}-\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r\frac{L_{\mathrm{<figure-callout\; id="0"\; label="m\; L\; r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">m</figure-callout>}}}{L_r}& 0\\ {\mathrm{\&omega;}}_r\frac{L_m}{L_r}& -\frac{L_m}{L_r{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}\\ {\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}\\ {\mathrm{\&omega;}}_r\end{array}\right]---\left(11\right)$

Formula (9) and formula (11) are carried out discretization, and the asynchronous machine system random process that is used for the estimation of depression of order EKF rotating speed has just become:

$\left\{\begin{array}{c}x\left(k\right)=f(x(k-1),u(k-1),w\left(k\right))={A^{\&prime;}}_{d}x(k-1)+{B^{\&prime;}}_{d}u(k-1)+w\left(k\right)\\ y\left(k\right)={C^{\&prime;}}_{d}x\left(k\right)+v\left(k\right)\end{array}\right.---\left(12\right)$

Wherein k is the sequence number of each state variable sequence,

x(k)＝[Ψ _αr(k)，Ψ _βr(k)，ω _r(k)] ^T，u(k)＝[i _αs(k)，i _βs(k)] ^T，

$y\left(k\right)=\left[\begin{array}{c}{u}_{\mathrm{\&alpha;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}\left(k\right)\\ u_{\mathrm{\&beta;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}\left(k\right)\end{array}\right],$

${A^{\&prime;}}_d=\left[\begin{array}{ccc}1-\frac{T_s}{{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r(k-1)T_s& 0\\ {\mathrm{\&omega;}}_r(k-1)T_s& 1-\frac{T_s}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\\ 0& 0& 1\end{array}\right],$

${B^{\&prime;}}_d=\left[\begin{array}{cc}\frac{L_m{T}_s}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\\ 0& \frac{L_m{T}_s}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}\\ 0& 0\end{array}\right],$

${C^{\&prime;}}_d=\left[\begin{array}{ccc}-\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r\left(k\right)\frac{L_{\mathrm{<figure-callout\; id="0"\; label="m\; L\; r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">m</figure-callout>}}}{L_r}& 0\\ {\mathrm{\&omega;}}_r\left(k\right)\frac{L_m}{L_r}& -\frac{L_m}{L_r{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\end{array}\right],$

W (k) and v (k) are respectively process noise and measure noise, T _sBe the sampling period;

Variable among formula (8)-Shi (12), meaning of parameters are all identical with formula (1);

So far, through obtained to supply to expand the asynchronous machine Mathematical Modeling that the Kalman filtering algorithm uses with up conversion;

Second step just obtained asynchronous machine depression of order EKF rotating speed algorithm for estimating with formula (12) substitution expansion Kalman filtering algorithm, and step is following:

1) assignment procedure noise covariance battle array Q and measurement noise covariance battle array R;

2) init state error covariance matrix P and state estimation

$\hat{x}\left(0\right)=E\left[x\left(0\right)\right],$ $P\left(0\right)=E{[x\left(0\right)-\hat{x}\left(0\right)][x\left(0\right)-\hat{x}\left(0\right)]}^T---\left(13\right)$

Wherein

is the state initial value, and P (0) is a state error covariance matrix initial value;

3) in each sampling period (k=1,2,3 ... ∞)

1. x (k) is done status predication:

${\hat{x}}^{-}\left(k\right)={A^{\&prime;}}_d\left(k\right)\hat{x}(k-1)+{B^{\&prime;}}_{d}u(k-1)---\left(14\right)$

2. error covariance matrix is done for the first time and is estimated:

P ^-(k)＝G(k)P(k-1)G(k) ^T+Q (15)

3. calculating K alman gain matrix K (k):

K(k)＝P-(k)M(k) ^T(M(k)P ^-(k)M(k) ^T+R) ^-1 (16)

4. correcting state is estimated

$\hat{x}\left(k\right)={\hat{x}}^{-}\left(k\right)+K\left(k\right)(y\left(k\right)-{C^{\&prime;}}_d{\hat{x}}^{-}\left(k\right))---\left(17\right)$

5. upgrade error covariance matrix P (k):

P(k)＝(I-K(k)M(k))P ^-(k) (18)

Wherein

$G\left(k\right)={\frac{\&PartialD;f}{\&PartialD;x}|}_{x=\hat{x}(k-1)}$

$=\left[\begin{array}{ccc}1-\frac{T_s}{{\mathrm{\&tau;}}_r}& -{\widehat{\mathrm{\&omega;}}}_r(k-1)T_s& -{\widehat{\mathrm{\&Psi;}}}_{\mathrm{\&beta;r}}(k-1)T_s\\ {\widehat{\mathrm{\&omega;}}}_r(k-1)T_s& 1-\frac{T_s}{{\mathrm{\&tau;}}_r}& {\widehat{\mathrm{\&Psi;}}}_{\mathrm{\&alpha;r}}(k-1)T_s\\ 0& 0& 1\end{array}\right]---\left(19\right)$

$M\left(k\right)={\frac{\&PartialD;h}{\&PartialD;x}|}_{x={\hat{x}}^{-}\left(k\right)}=\left[\begin{array}{ccc}-\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& -{{\widehat{\mathrm{\&omega;}}}^{-}}_r\left(k\right)\frac{L_m}{L_r}& -{{\widehat{\mathrm{\&Psi;}}}^{-}}_{\mathrm{\&beta;r}}\left(k\right)\frac{L_m}{L_r}\\ {{\widehat{\mathrm{\&omega;}}}^{-}}_r\left(k\right)\frac{L_m}{L_r}& -\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& {{\widehat{\mathrm{\&Psi;}}}^{-}}_{\mathrm{\&alpha;r}}\left(k\right)\frac{L_m}{L_r}\end{array}\right]---\left(20\right)$

Subscript "-" among formula (13)-Shi (20) is all represented premeasuring, and subscript " ∧ " is all represented estimator; Variable, meaning of parameters are all identical with formula (1);

So far, obtained asynchronous machine depression of order EKF rotating speed algorithm for estimating at each sampling period T _sInterior concrete calculation procedure and computing formula are for the hardware language description of next step FPGA is got ready;

In the 3rd step, design FPGA realizes the algorithm structure that depression of order EKF rotating speed is estimated:

According to formula (14), control u (k) multiply by B ' _dBattle array 96 backs and state transitions battle array A ' _d(k) and state Long-pending 97 obtain the result that 1. EKF the goes on foot status predication through adder 95 additions Wherein, u (k) is the stator current vector; Wherein

It is the state estimation that a last bat obtains

Obtain through 98 delays, one bat of unit delay unit; Predicted state

With a measurement battle array C ' _d(k) 91 multiply each other, subtracter 92 be used for calculating measuring amount y (k) therewith the difference of product give Kalman gain adjustment module 93, i.e. the of EKF 3. step; Before this, FPGA should trigger the 1. the step calculate in the computing circuit of trigger error Estimates on Covariance Matrix, upgrade P by formula (15) ^-(k), promptly accomplish the 2. step of EKF; Use 93 pairs of predicted states of Kalman gain adjustment module

Proofread and correct by formula 17, adder 94 is used on predicted state, adding the correcting value by the long-pending decision of Kalman gain and measure error, obtains final state estimation

Promptly accomplish the 4. step of EKF; After this, the computing circuit that FPGA answers the trigger error covariance matrix to upgrade upgrades P (k) by formula (18), promptly accomplishes the 5. step of EKF; So far, all interior calculating of sampling period of EKF have been accomplished;

The computing formula of measuring amount y (k) is following:

$y\left(k\right)=\left[\begin{array}{c}{u}_{\mathrm{\&alpha;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}\left(k\right)\\ u_{\mathrm{\&beta;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}\left(k\right)\end{array}\right]---\left(21\right)$

Plus and minus calculation in subtracter 991 and subtracter 992 perfects (21), stator current measured value i _S(k)=[i _{α s}(k), i _{β s}(k)] ^TMultiply each other obtain long-pending and the differential signal of the stator current that calculates by weighted network 994, coefficient 996 with coefficient 993

With the product of coefficient 995 as subtrahend and stator voltage command value

Poor, measuring amount y (k) just; Differential signal to stator current

Calculating adopt method---the three point value differential of evaluation differential in the numerical analysis:

${\stackrel{\&CenterDot;}{i}}_s\left(k\right)=\frac{1}{2T_s}(i_s(k-2)-4i_s(k-1)+3i_s\left(k\right))---\left(22\right)$

Displacement weighted sum in the following formula (22) is accomplished i by weighted network 994 _S(k) each is clapped and is postponed to realize that by a unit delay unit this weighted sum and coefficient 996 multiply each other and be the differential of stator current

The 4th step,

Based on the algorithm structure of above description, FPGA is carried out hardware language VHDL describes, the equal using system clock of all modules CLK is as clock signal, Rstx signal that each module receives and Cwel signal deciding their calculating sequential; The Rstx signal is low effective enable signal, and Cwel is a register type, count range 0 to 12, and the duration of each count value is 3 CLK, the matrix multiplication that is designed can be accomplished at 3 CLK; In each count cycle, when Rst1 when low, measuring amount computing module 937, status predication module 932 and error covariance matrix estimation module 936 are enabled, computing begins, the operation time that error covariance matrix estimation module 936 needs is the longest, is 18 CLK; After 18 CLK; Be that the Cwel count value is after 6; 1. the of the calculating measuring amount computing module 937 of measuring amount y (k) and EKF go on foot status predication module 932 and the and 2. go on foot error covariance matrix estimation module 936 and all carried out computing concomitantly, for the of EKF 3. go on foot Kalman gain calculating 935 all set when time each variate-value; The Cwel count value is 7～9 these 3 count values, and Rst3 is when low, and Kalman gain calculating 935 is enabled; The Cwel count value is 10～12 these 3 count values; And when Rst4 is low; State correction module 933 and error covariance matrix update module 934 are enabled simultaneously; The Cwel count value is after 12, and 4. of EKF goes on foot state correction module 933 and and 5. go on foot error covariance matrix update module 934 and all carried out computing concomitantly, as inferior state estimation value i _{α s}, i _{β s}, Ψ _{α r}, Ψ _{β r}, ω _rBe ready to be sent to the parallel port and passed master control DSP back; i _{α s}, i _{β s}Be the estimated value of α, β axle stator current, Ψ _{α r}, Ψ _{β r}Be the estimated value of α, β axle rotor flux, ω _rEstimated value for rotor electrical angle speed.

The invention has the beneficial effects as follows; The FPGA that utilization has operation core realizes EKF rotating speed algorithm for estimating concurrently with electronic circuit; This can bring the benefit of two aspects: 1) greatly reduce master control DSP the real-time operation amount on pressure; For speed, Current Control stay more memory space and computational space, make more senior speed ring, electric current loop algorithm on DSP, realize becoming possibility; 2) can in the time of 1 μ s level, accomplish with the EKF rotating speed algorithm for estimating of FPGA Parallel Implementation; Therefore the sampling period of EKF algorithm mainly by and master control DSP between the data transmission bauds decision; Make the rotating speed algorithm for estimating can select the less sampling period, thereby the rotating speed estimated accuracy has improved greatly.

Description of drawings

Fig. 1 is that the system of the direct vector control apparatus of realization asynchronous machine Speedless sensor of existing routine realizes block diagram;

Fig. 2 is a system block diagram of the present invention, and the depression of order EKF algorithm 90 that wherein is used for real-time rotating speed estimation is realized with the form of hardwareization by FPGA 130;

Fig. 3 is the algorithm structure sketch map that the depression of order EKF rotating speed that adopts among the present invention is estimated;

Fig. 4 is the hardware description top layer project organization of depression of order EKF rotating speed algorithm for estimating FPGA Parallel Implementation.

Embodiment

Fig. 1 is that the system of the direct vector control apparatus of realization asynchronous machine Speedless sensor of existing routine realizes block diagram; Wherein contain EKF rotating speed algorithm for estimating 90 and all use dsp programs to realize in interior whole vector control algorithm 120, promptly among the figure whole vector control algorithm 120 need take DSP the memory space resource and running time resource.The use EKF of this routine does the rotating speed estimation provides the direct vector control apparatus of asynchronous machine of speed feedback to comprise the closed-loop control of speed, torque current and rotor flux exciting current.Speed pi regulator 30 provides torque reference value according to the deviation between rotary speed instruction

and the speed feedback value

; Under the constant situation of rotor flux; Torque is proportional to the stator current torque component; Therefore; The output of speed regulator 30 provides with the form of torque current reference value

; Make the rotating speed

of the asynchronous machine 10 that estimates from extended Kalman filter 90 can follow the tracks of rotary speed instruction value

deviation signal

and

are input to current regulator 40 separately; Diaxon voltage reference value

and

extended Kalman filter 90 of obtaining in the field orientation dq axle system provide unit vector 110 for vector control system; The space phase

that indicates the rotor flux vector that estimates through this unit vector is directed coordinate with this phase place as the dq axle; Make dq axle system directed along the rotor flux direction; Here it is said rotor field-oriented control, i.e. vector control.At this moment, program 60 is sent 6 road pwm signals driving three-phase voltage-type inverter 70 according to these two voltage reference value components switching device takes place through voltage reference value

and

SVPWM that coordinate transform 50 is transformed to voltage reference value on the rest frame α β in the diaxon voltage reference value

in the dq axle system and

.Current Hall detector 20 is used for detecting the stator current i of asynchronous machine 10 _As, i _Bs, i _Cs, after 80s through three-phase/two phase inversion, obtain the stator current detected value i on the rest frame α β _{α s}And i _{β s}EKF 90 utilizes the stator voltage command value

With

And stator current detected value i _{α s}And i _{β s}For Direct Vector Control provides the estimated value of stator current, rotor flux and unit vector, wherein stator current estimated value With

Be transformed to two shaft current estimated values in the field orientation dq axle system through coordinate transform 100

With

Value of feedback and torque current reference value as current inner loop 40

And exciting current reference value

Relatively produce deviation respectively as the input of current regulator 43,44.

The shortcoming of the conventional direct vector control apparatus of asynchronous machine Speedless sensor is vector control, rotating speed algorithm for estimating all with the dsp program realization, and this speed of service to DSP has proposed very high requirement.Through actual measurement, using the running time of EKF rotating speed algorithm for estimating on TMS320F2812 fixed DSP (dominant frequency 150MHz) of 32 multiplication is 50 μ s, and this makes the computing of in rational control cycle, accomplishing the current inner loop control algolithm hardly maybe.On the other hand; The arithmetic speed of DSP is difficult to support realize the high-accuracy speed estimation with EKF at present; The estimated accuracy of EKF improves rapidly with the shortening in sampling period; Want to be met the rotating speed estimated accuracy of commercial Application index, 50 μ s are the maximum critical value in EKF algorithm sampling period basically.In fact, the suitable EKF algorithm sampling period should be about 10 μ s.But this instruction cycle that just requires DSP is more than 500MHz.Why present this also is still mainly to adopt one of reason of model reference adaptive rotating speed algorithm for estimating in the vector control type frequency converter.But as before carried, the model reference adaptive method for estimating rotating speed receives the influence of initial value for integral and drift very big, the rotating speed estimated result is inaccurate, and poor-performing during low speed.

The present invention provides the pure electronic type speed feedback method in the asynchronous machine speed-less sensor vector control system; Consider the level of hardware of present FPGA, the present invention temporarily uses the EKF algorithm of depression of order to illustrate as illustrative examples how to use this method that the asynchronous machine rotating speed is estimated in real time.System block diagram of the present invention is as shown in Figure 2; Be that with the difference of the asynchronous machine Speedless sensor Direct Vector Control of routine shown in Figure 1 the EKF rotating speed of being realized by the dsp software programming is originally estimated that 90 parts change by the FPGA hardware circuit to be realized, has alleviated the computation burden of DSP.For improving the data communication speed between DSP and the FPGA, use parallel port 140 and parallel port 150 transmission data, if the parallel port resource of DSP is less, then can use SPI to carry out Serial Data Transfer Mode.The estimation 90 of depression of order EKF rotating speed is estimated the rotating speed of asynchronous machine according to algorithm according to stator voltage command value and stator current measured value that DSP provides in real time; And the stator current of estimating, rotor flux and rotating speed passed to master control DSP again back the reading in, write out of data through the parallel port 150 of FPGA by 160 controls of the parallel port data control logic among the FPGA.The algorithm structure of depression of order EKF rotating speed estimation 90 is seen Fig. 3.Fig. 4 is seen in the FPGA top layer design of depression of order EKF rotating speed estimation 90.

Concrete feedback method is following:

The first step obtains the discrete time formal representation of asynchronous machine Mathematical Modeling:

Here be estimated as illustrative examples with the asynchronous machine EKF rotating speed behind the depression of order, at first according to two mutually static α β shaft models of asynchronous machine:

$\left[\begin{array}{c}{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}\\ {\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&alpha;r}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&beta;r}}\end{array}\right]=\left[\begin{array}{cccc}-(\frac{R_s}{\mathrm{\&sigma;}{L}_s}+\frac{1-\mathrm{\&sigma;}}{\mathrm{\&sigma;}{\mathrm{\&tau;}}_r})& 0& \frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r{\mathrm{\&tau;}}_r}& {\mathrm{\&omega;}}_r\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r}\\ 0& -(\frac{R_s}{\mathrm{\&sigma;}{L}_s}+\frac{1-\mathrm{\&sigma;}}{\mathrm{\&sigma;}{\mathrm{\&tau;}}_r})& -{\mathrm{\&omega;}}_r\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r}& \frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r{\mathrm{\&tau;}}_r}\\ \frac{L_m}{{\mathrm{\&tau;}}_r}& 0& -\frac{1}{{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r\\ 0& \frac{L_m}{{\mathrm{\&tau;}}_r}& {\mathrm{\&omega;}}_r& -\frac{1}{{\mathrm{\&tau;}}_r}\end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{\&alpha;s}}\\ i_{\mathrm{\&beta;s}}\\ {\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}\\ {\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}\end{array}\right]+$

$\left[\begin{array}{cc}\frac{1}{\mathrm{\&sigma;}{L}_s}& 0\\ 0& \frac{1}{\mathrm{\&sigma;}{L}_s}\\ 0& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{u}_{\mathrm{\&alpha;s}}\\ u_{\mathrm{\&beta;s}}\end{array}\right]---\left(1\right)$

In the formula (1), R _s, R _rBe stator and rotor resistance parameters; L _s, L _rBe the stator and rotor inductance; L _mBe mutual inductance; σ is a magnetic leakage factor, σ=1-(L _m ²/ L _sL _r); τ _rBe rotor time constant, τ _r=L _r/ R _rω _rBe rotor electrical angle speed; u _{α s}, u _{β s}Be α, β axle stator voltage; i _{α s}, i _{β s}Be α, β axle stator current; Ψ _{α r}, Ψ _{β r}Be α, β axle rotor flux.Subscript ". " expression is differentiated, like

Differential equation group form with formula (1) writing formula (2)～(5):

${\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}=-(\frac{R_s}{\mathrm{\&sigma;}{L}_s}+\frac{1-\mathrm{\&sigma;}}{\mathrm{\&sigma;}{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}+\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}+{\mathrm{\&omega;}}_r\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r}{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}+\frac{1}{\mathrm{\&sigma;}{L}_s}{u}_{\mathrm{\&alpha;s}}---\left(2\right)$

${\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}=-(\frac{R_s}{\mathrm{\&sigma;}{L}_s}+\frac{1-\mathrm{\&sigma;}}{\mathrm{\&sigma;}{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}-{\mathrm{\&omega;}}_r\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}+\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}+\frac{1}{\mathrm{\&sigma;}{L}_s}{u}_{\mathrm{\&beta;s}}---\left(3\right)$

${\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&alpha;r}}=\frac{L_m}{{\mathrm{\&tau;}}_r}{i}_{\mathrm{\&alpha;s}}-\frac{1}{{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}-{\mathrm{\&omega;}}_r{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}---\left(4\right)$

${\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&beta;r}}=\frac{L_m}{{\mathrm{\&tau;}}_r}{i}_{\mathrm{\&beta;s}}+{\mathrm{\&omega;}}_r{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}-\frac{1}{{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}---\left(5\right)$

Formula (2), (3) both sides are with taking advantage of σ L _s,

$\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}=-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}+\frac{L_m}{L_r{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}+{\mathrm{\&omega;}}_r\frac{L_m}{L_r}{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}+u_{\mathrm{\&alpha;s}}---\left(6\right)$

$\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}=-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}-{\mathrm{\&omega;}}_r\frac{L_m}{L_r}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}+\frac{L_m}{L_r{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}+u_{\mathrm{\&beta;s}}---\left(7\right)$

Suppose that change in rotational speed is very little in a sampling period, can use approx

Replace mechanical equation, reelect x=[Ψ _{α r}, Ψ _{β r}, ω _r] ^TMake state, u=[i _{α s}, i _{β s}] ^TControl, the transposition computing of subscript T representing matrix is deformed into output equation with formula (6), (7) on this basis.After the arrangement, the 3 scalariform attitude equations that obtain the asynchronous machine system are:

$\stackrel{\&CenterDot;}{x}=f(x,u)---\left(8\right)$

Promptly

$\left[\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&alpha;r}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&beta;r}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_r\end{array}\right]=\left[\begin{array}{ccc}-\frac{1}{{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}& 0\\ {\mathrm{\&omega;}}_r& -\frac{1}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}\\ {\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}\\ {\mathrm{\&omega;}}_r\end{array}\right]+\left[\begin{array}{cc}\frac{L_m}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\\ 0& \frac{L_m}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}\\ 0& 0\end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{\&alpha;s}}\\ i_{\mathrm{\&beta;s}}\end{array}\right]---\left(9\right)$

The measurement equation is:

y＝h(x) (10)

Promptly

$y=\left[\begin{array}{c}{u}_{\mathrm{\&alpha;s}}-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}\\ u_{\mathrm{\&beta;s}}-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}\end{array}\right]=\left[\begin{array}{ccc}-\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r\frac{L_{\mathrm{<figure-callout\; id="0"\; label="m\; L\; r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">m</figure-callout>}}}{L_r}& 0\\ {\mathrm{\&omega;}}_r\frac{L_m}{L_r}& -\frac{L_m}{L_r{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}\\ {\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}\\ {\mathrm{\&omega;}}_r\end{array}\right]---\left(11\right)$

Discretization is carried out in (9) and (11), and the asynchronous machine system random process that is used for the estimation of depression of order EKF rotating speed has just become:

$\left\{\begin{array}{c}x\left(k\right)=f(x(k-1),u(k-1),w\left(k\right))={A^{\&prime;}}_{d}x(k-1)+{B^{\&prime;}}_{d}u(k-1)+w\left(k\right)\\ y\left(k\right)={C^{\&prime;}}_{d}x\left(k\right)+v\left(k\right)\end{array}\right.---\left(12\right)$

Wherein k is the sequence number of each variable sequence,

x(k)＝[Ψ _αr(k)，Ψ _βr(k)，ω _r(k)] ^T，u(k)＝[i _αs(k)，i _βs(k)] ^T，

$y\left(k\right)=\left[\begin{array}{c}{u}_{\mathrm{\&alpha;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}\left(k\right)\\ u_{\mathrm{\&beta;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}\left(k\right)\end{array}\right],$

${A^{\&prime;}}_d=\left[\begin{array}{ccc}1-\frac{T_s}{{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r(k-1)T_s& 0\\ {\mathrm{\&omega;}}_r(k-1)T_s& 1-\frac{T_s}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\\ 0& 0& 1\end{array}\right],$

${B^{\&prime;}}_d=\left[\begin{array}{cc}\frac{L_m{T}_s}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\\ 0& \frac{L_m{T}_s}{{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}\\ 0& 0\end{array}\right],$

${C^{\&prime;}}_d=\left[\begin{array}{ccc}-\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r\left(k\right)\frac{L_{\mathrm{<figure-callout\; id="0"\; label="m\; L\; r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">m</figure-callout>}}}{L_r}& 0\\ {\mathrm{\&omega;}}_r\left(k\right)\frac{L_m}{L_r}& -\frac{L_m}{L_r{\mathrm{\&tau;}}_{\mathrm{<figure-callout\; id="0"\; label="r"\; filenames="HDA0000148136900000021.png"\; state="\{\{state\}\}">r</figure-callout>}}}& 0\end{array}\right],$ W (k) and v (k) are respectively process noise and measure noise, T _sBe the sampling period.

Variable among formula (8)-Shi (12), meaning of parameters are all identical with formula (1).

So far, through obtained to supply to expand the asynchronous machine Mathematical Modeling that the Kalman filtering algorithm uses with up conversion.

Second step just obtained asynchronous machine depression of order EKF rotating speed algorithm for estimating with formula (12) substitution expansion Kalman filtering algorithm, and step is following:

1) assignment procedure noise covariance battle array Q and measurement noise covariance battle array R;

2) init state error covariance matrix P and state estimation

$\hat{x}\left(0\right)=E\left[x\left(0\right)\right],$ $P\left(0\right)=E{[x\left(0\right)-\hat{x}\left(0\right)][x\left(0\right)-\hat{x}\left(0\right)]}^T---\left(13\right)$

Wherein is the state initial value, and P (0) is a state error covariance matrix initial value.

3) in each sampling period (k=1,2,3 ... ∞)

1. x (k) is done status predication:

${\hat{x}}^{-}\left(k\right)={A^{\&prime;}}_d\left(k\right)\hat{x}(k-1)+{B^{\&prime;}}_{d}u(k-1)---\left(14\right)$

2. error covariance matrix is done for the first time and is estimated:

P ^-(k)＝G(k)P(k-1)G(k) ^T+Q (15)

3. calculating K alman gain matrix K (k):

K(k)＝P ^-(k)M(k) ^T(M(k)P ^-(k)M(k) ^T+R) ^-1 (16)

4. correcting state is estimated

$\hat{x}\left(k\right)={\hat{x}}^{-}\left(k\right)+K\left(k\right)(y\left(k\right)-{C^{\&prime;}}_d{\hat{x}}^{-}\left(k\right))---\left(17\right)$

5. upgrade error covariance matrix P (k):

P(k)＝(I-K(k)M(k))P ^-(k) (18)

Wherein

$G\left(k\right)={\frac{\&PartialD;f}{\&PartialD;x}|}_{x=\hat{x}(k-1)}$

$=\left[\begin{array}{ccc}1-\frac{T_s}{{\mathrm{\&tau;}}_r}& -{\widehat{\mathrm{\&omega;}}}_r(k-1)T_s& -{\widehat{\mathrm{\&Psi;}}}_{\mathrm{\&beta;r}}(k-1)T_s\\ {\widehat{\mathrm{\&omega;}}}_r(k-1)T_s& 1-\frac{T_s}{{\mathrm{\&tau;}}_r}& {\widehat{\mathrm{\&Psi;}}}_{\mathrm{\&alpha;r}}(k-1)T_s\\ 0& 0& 1\end{array}\right]---\left(19\right)$

$M\left(k\right)={\frac{\&PartialD;h}{\&PartialD;x}|}_{x={\hat{x}}^{-}\left(k\right)}=\left[\begin{array}{ccc}-\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& -{{\widehat{\mathrm{\&omega;}}}^{-}}_r\left(k\right)\frac{L_m}{L_r}& -{{\widehat{\mathrm{\&Psi;}}}^{-}}_{\mathrm{\&beta;r}}\left(k\right)\frac{L_m}{L_r}\\ {{\widehat{\mathrm{\&omega;}}}^{-}}_r\left(k\right)\frac{L_m}{L_r}& -\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& {{\widehat{\mathrm{\&Psi;}}}^{-}}_{\mathrm{\&alpha;r}}\left(k\right)\frac{L_m}{L_r}\end{array}\right]---\left(20\right)$

Subscript "-" among formula (13)-Shi (20) is all represented premeasuring, and subscript " ∧ " is all represented estimator.Variable, meaning of parameters are all identical with formula (1).

So far, obtained asynchronous machine depression of order EKF rotating speed algorithm for estimating at each sampling period T _sInterior concrete calculation procedure and computing formula are for the hardware language description of next step FPGA is got ready.

In the 3rd step, design FPGA realizes the algorithm structure that depression of order EKF rotating speed is estimated:

During the algorithm described with following formula (14)～(20) with the FPGA practical implementation, algorithm structure such as Fig. 3.According to (14) formula, control u (k) (being the stator current vector here) multiply by B ' _dBattle array 96 backs and state transitions battle array A ' _d(k) and state

Long-pending 97 obtain the result that 1. EKF the goes on foot status predication through adder 95 additions

Wherein

It is the state estimation that a last bat obtains

Obtain through 98 delays, one bat of unit delay unit.Predicted state

With a measurement battle array C ' _d(k) 91 multiply each other, subtracter 92 be used for calculating measuring amount y (k) therewith the difference of product give Kalman gain adjustment module 93, i.e. the of EKF 3. step.Before this, FPGA should trigger the 1. the step calculate in the computing circuit of trigger error Estimates on Covariance Matrix, upgrade P by formula (15) ^-(k), promptly accomplish the 2. step of EKF.Use 93 pairs of predicted states of Kalman gain adjustment module to proofread and correct by formula (17); Adder 94 is used on predicted state adding the correcting value by the long-pending decision of Kalman gain and measure error, obtains final state estimation

and promptly accomplishes the of EKF and 4. go on foot.After this, the computing circuit that FPGA answers the trigger error covariance matrix to upgrade upgrades P (k) by formula (18), promptly accomplishes the 5. step of EKF.So far, all interior calculating of sampling period of EKF have been accomplished.

Need to prove, when the measuring amount y (k) that the present invention carries out choosing when the EKF rotating speed is estimated need begin in each sampling period of EKF according to the stator voltage command value

And stator current measured value i _S(k)=[i _{α s}(k), i _{β s}(k)] ^TCalculate and obtain.Can be found out that by formula (11) measuring y no longer is simple state variable, the observation battle array also no longer is linear.Calculating 99 algorithm structures to this part calculating below in conjunction with the measuring amount of Fig. 3 specifies.The computing formula of measuring amount y (k) is following:

$y\left(k\right)=\left[\begin{array}{c}{u}_{\mathrm{\&alpha;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}\left(k\right)\\ u_{\mathrm{\&beta;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}\left(k\right)\end{array}\right]---\left(21\right)$

Plus and minus calculation in subtracter 991 and subtracter 992 perfects (21), stator current measured value i _S(k) multiply each other obtain long-pending with coefficient 993 and the differential signal of the stator current that calculates by weighted network 994, coefficient 996

With the product of coefficient 995 as subtrahend and stator voltage command value

Poor, be exactly the measuring amount y (k) that uses among the present invention.The calculating of the differential signal

of stator current is adopted method---the three point value differential of evaluation differential in the numerical analysis:

${\stackrel{\&CenterDot;}{i}}_s\left(k\right)=\frac{1}{2T_s}(i_s(k-2)-4i_s(k-1)+3i_s\left(k\right))---\left(22\right)$

Displacement weighted sum in the following formula is accomplished i by weighted network 994 _S(k) each is clapped and is postponed to realize that by a unit delay unit this weighted sum and coefficient 996 multiply each other and be the differential of stator current

The 4th step, based on the algorithm structure of above description, FPGA to be carried out hardware language VHDL describe (promptly in the FPGA of PC development environment, programming), the top layer design is as shown in Figure 4.The equal using system clock of all modules CLK is as clock signal, Rstx signal that each module receives and Cwel signal deciding their calculating sequential.The Rstx signal is low effective enable signal, and Cwel is a register type, count range 0 to 12, and the duration of each count value is 3 CLK, the matrix multiplication that is designed can be accomplished at 3 CLK.In each count cycle (being the EKF sampling period); When Rst1 when low, measuring amount computing module 937, status predication module 932 and error covariance matrix estimation module 936 are enabled, computing begins; The operation time that error covariance matrix estimation module 936 needs is the longest, is 18 CLK.After 18 CLK, promptly the Cwel count value is after 6, and 1. the of the calculating 937 of measuring amount y (k) and EKF go on foot 932 and the and 2. go on foot 936 and all carried out computing concomitantly, for the of EKF 3. go on foot Kalman gain calculating 935 all set when time each variate-value.The Cwel count value is 7～9 these 3 count values, and Rst3 is when low, and Kalman gain calculating 935 is enabled.The Cwel count value is 10～12 these 3 count values; And when Rst4 was low, state correction module 933 and error covariance matrix update module 934 were enabled simultaneously, and the Cwel count value is after 12; 4. of EKF goes on foot 933 and and 5. goes on foot 934 and all carried out computing concomitantly, as inferior state estimation value i _{α s}, i _{β s}, Ψ _{α r}, Ψ _{β r}, ω _rBe ready to be sent to the parallel port and passed master control DSP back; i _{α s}, i _{β s}Be the estimated value of α, β axle stator current, Ψ _{α r}, Ψ _{β r}Be the estimated value of α, β axle rotor flux, ω _rEstimated value for rotor electrical angle speed.

Use the change of FPGA hardware time of implementation in asynchronous machine depression of order EKF rotating speed each sampling period of algorithm for estimating be 36 CLK, be example with the EP3C16Q240C8 of Altera, when the external crystal-controlled oscillation frequency was 30MHz, the time of implementation of algorithm was merely 1.2 μ s.So short algorithm execution time makes selects the EKF rotating speed algorithm for estimating sampling period of the 10 μ s orders of magnitude to become possibility, has therefore improved estimated accuracy greatly.

Concerning master control DSP; A large amount of complex calculation that have degree of precision to require are independently born by FPGA; No longer with DSP in other complicated control algolithm competes for memory resource of high-performance and real time execution time resource, make the more high performance speed of asynchronous machine, torque control become possibility.

The present invention realizes the method for the online estimation of asynchronous machine rotating speed with pure hardware circuit; The main body of this hardware circuit is FPGA (Field Programmable Gate Array-field programmable gate array) chip; Calculation function in this chip is defined as expansion Kalman (Kalman) the filtering turn count algorithm of depression of order by the developer with the form of hardware description language especially, so just can in the asynchronous machine speed-less sensor vector control system, be independent of master control DSP (Digital Signal Processor-digital signal processor) chip and provide asynchronous machine high accuracy speed feedback information in real time.

Claims (1)

1. the pure electronic type speed feedback of asynchronous machine method is characterized in that, the concrete operations step is following:

The first step obtains the discrete time formal representation of asynchronous machine Mathematical Modeling:

At first according to two mutually static α β shaft models of asynchronous machine:

$\left[\begin{array}{c}{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}\\ {\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&alpha;r}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&beta;r}}\end{array}\right]=\left[\begin{array}{cccc}-(\frac{R_s}{\mathrm{\&sigma;}{L}_s}+\frac{1-\mathrm{\&sigma;}}{\mathrm{\&sigma;}{\mathrm{\&tau;}}_r})& 0& \frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r{\mathrm{\&tau;}}_r}& {\mathrm{\&omega;}}_r\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r}\\ 0& -(\frac{R_s}{\mathrm{\&sigma;}{L}_s}+\frac{1-\mathrm{\&sigma;}}{\mathrm{\&sigma;}{\mathrm{\&tau;}}_r})& -{\mathrm{\&omega;}}_r\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r}& \frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r{\mathrm{\&tau;}}_r}\\ \frac{L_m}{{\mathrm{\&tau;}}_r}& 0& -\frac{1}{{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r\\ 0& \frac{L_m}{{\mathrm{\&tau;}}_r}& {\mathrm{\&omega;}}_r& -\frac{1}{{\mathrm{\&tau;}}_r}\end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{\&alpha;s}}\\ i_{\mathrm{\&beta;s}}\\ {\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}\\ {\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}\end{array}\right]+$

$\left[\begin{array}{cc}\frac{1}{\mathrm{\&sigma;}{L}_s}& 0\\ 0& \frac{1}{\mathrm{\&sigma;}{L}_s}\\ 0& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{u}_{\mathrm{\&alpha;s}}\\ u_{\mathrm{\&beta;s}}\end{array}\right]---\left(1\right)$

In the formula (1), R _s, R _rBe stator and rotor resistance parameters; L _s, L _rBe the stator and rotor inductance; L _mBe mutual inductance; σ is a magnetic leakage factor, σ=1-(L _m ²/ L _sL _r); τ _rBe rotor time constant, τ _r=L _r/ R _rω _rBe rotor electrical angle speed; u _{α s}, u _{β s}Be α, β axle stator voltage; i _{α s}, i _{β s}Be α, β axle stator current; Ψ _{α r}, Ψ _{β r}Be α, β axle rotor flux;

Differential equation group form with formula (1) writing formula (2)～(5):

${\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}=-(\frac{R_s}{\mathrm{\&sigma;}{L}_s}+\frac{1-\mathrm{\&sigma;}}{\mathrm{\&sigma;}{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}+\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}+{\mathrm{\&omega;}}_r\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r}{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}+\frac{1}{\mathrm{\&sigma;}{L}_s}{u}_{\mathrm{\&alpha;s}}---\left(2\right)$

${\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}=-(\frac{R_s}{\mathrm{\&sigma;}{L}_s}+\frac{1-\mathrm{\&sigma;}}{\mathrm{\&sigma;}{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}-{\mathrm{\&omega;}}_r\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}+\frac{L_m}{\mathrm{\&sigma;}{L}_s{L}_r{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}+\frac{1}{\mathrm{\&sigma;}{L}_s}{u}_{\mathrm{\&beta;s}}---\left(3\right)$

${\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&alpha;r}}=\frac{L_m}{{\mathrm{\&tau;}}_r}{i}_{\mathrm{\&alpha;s}}-\frac{1}{{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}-{\mathrm{\&omega;}}_r{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}---\left(4\right)$

${\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&beta;r}}=\frac{L_m}{{\mathrm{\&tau;}}_r}{i}_{\mathrm{\&beta;s}}+{\mathrm{\&omega;}}_r{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}-\frac{1}{{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}---\left(5\right)$

Formula (2), (3) both sides are with taking advantage of σ L _s,

$\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}=-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}+\frac{L_m}{L_r{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}+{\mathrm{\&omega;}}_r\frac{L_m}{L_r}{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}+u_{\mathrm{\&alpha;s}}---\left(6\right)$

$\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}=-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}-{\mathrm{\&omega;}}_r\frac{L_m}{L_r}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}+\frac{L_m}{L_r{\mathrm{\&tau;}}_r}{\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}+u_{\mathrm{\&beta;s}}---\left(7\right)$

Suppose that change in rotational speed is very little in a sampling period, can use approx

Replace mechanical equation, reelect x=[Ψ _{α r}, Ψ _{β r}, ω _r] ^TMake state, u=[i _{α s}, i _{β s}] ^TControl, the transposition computing of subscript T representing matrix is deformed into output equation with formula (6), formula (7) on this basis; After the arrangement, the 3 scalariform attitude equations that obtain the asynchronous machine system are:

$\stackrel{\&CenterDot;}{x}=f(x,u)---\left(8\right)$

Promptly

$\left[\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&alpha;r}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&Psi;}}}_{\mathrm{\&beta;r}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_r\end{array}\right]=\left[\begin{array}{ccc}-\frac{1}{{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r& 0\\ {\mathrm{\&omega;}}_r& -\frac{1}{{\mathrm{\&tau;}}_r}& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}\\ {\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}\\ {\mathrm{\&omega;}}_r\end{array}\right]+\left[\begin{array}{cc}\frac{L_m}{{\mathrm{\&tau;}}_r}& 0\\ 0& \frac{L_m}{{\mathrm{\&tau;}}_r}\\ 0& 0\end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{\&alpha;s}}\\ i_{\mathrm{\&beta;s}}\end{array}\right]---\left(9\right)$

The measurement equation is:

y＝h(x) (10)

Promptly

$y=\left[\begin{array}{c}{u}_{\mathrm{\&alpha;s}}-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}\\ u_{\mathrm{\&beta;s}}-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}\end{array}\right]=\left[\begin{array}{ccc}-\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r\frac{L_m}{L_r}& 0\\ {\mathrm{\&omega;}}_r\frac{L_m}{L_r}& -\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Psi;}}_{\mathrm{\&alpha;r}}\\ {\mathrm{\&Psi;}}_{\mathrm{\&beta;r}}\\ {\mathrm{\&omega;}}_r\end{array}\right]---\left(11\right)$

Formula (9) and formula (11) are carried out discretization, and the asynchronous machine system random process that is used for the estimation of depression of order EKF rotating speed has just become:

$\left\{\begin{array}{c}x\left(k\right)=f(x(k-1),u(k-1),w\left(k\right))={A^{\&prime;}}_{d}x(k-1)+{B^{\&prime;}}_{d}u(k-1)+w\left(k\right)\\ y\left(k\right)={C^{\&prime;}}_{d}x\left(k\right)+v\left(k\right)\end{array}\right.---\left(12\right)$

Wherein k is the sequence number of each state variable sequence,

x(k)＝[Ψ _αr(k)，Ψ _βr(k)，ω _r(k)] ^T，u(k)＝[i _αs(k)，i _βs(k)] ^T，

$y\left(k\right)=\left[\begin{array}{c}{u}_{\mathrm{\&alpha;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}\left(k\right)\\ u_{\mathrm{\&beta;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}\left(k\right)\end{array}\right],$

${A^{\&prime;}}_d=\left[\begin{array}{ccc}1-\frac{T_s}{{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r(k-1)T_s& 0\\ {\mathrm{\&omega;}}_r(k-1)T_s& 1-\frac{T_s}{{\mathrm{\&tau;}}_r}& 0\\ 0& 0& 1\end{array}\right],$

${B^{\&prime;}}_d=\left[\begin{array}{cc}\frac{L_m{T}_s}{{\mathrm{\&tau;}}_r}& 0\\ 0& \frac{L_m{T}_s}{{\mathrm{\&tau;}}_r}\\ 0& 0\end{array}\right],$

${C^{\&prime;}}_d=\left[\begin{array}{ccc}-\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& -{\mathrm{\&omega;}}_r\left(k\right)\frac{L_m}{L_r}& 0\\ {\mathrm{\&omega;}}_r\left(k\right)\frac{L_m}{L_r}& -\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& 0\end{array}\right],$

W (k) and v (k) are respectively process noise and measure noise, T _sBe the sampling period;

Variable among formula (8)-Shi (12), meaning of parameters are all identical with formula (1);

So far, through obtained to supply to expand the asynchronous machine Mathematical Modeling that the Kalman filtering algorithm uses with up conversion;

Second step just obtained asynchronous machine depression of order EKF rotating speed algorithm for estimating with formula (12) substitution expansion Kalman filtering algorithm, and step is following:

1) assignment procedure noise covariance battle array Q and measurement noise covariance battle array R;

2) init state error covariance matrix P and state estimation

$\hat{x}\left(0\right)=E\left[x\left(0\right)\right],$ $P\left(0\right)=E{[x\left(0\right)-\hat{x}\left(0\right)][x\left(0\right)-\hat{x}\left(0\right)]}^T---\left(13\right)$

Wherein

is the state initial value, and P (0) is a state error covariance matrix initial value;

3) in each sampling period (k=1,2,3 ... ∞)

1. x (k) is done status predication:

${\hat{x}}^{-}\left(k\right)={A^{\&prime;}}_d\left(k\right)\hat{x}(k-1)+{B^{\&prime;}}_{d}u(k-1)---\left(14\right)$

2. error covariance matrix is done for the first time and is estimated:

P-(k)＝G(k)P(k-1)G(k) ^T+Q (15)

3. calculating K alman gain matrix K (k):

K(k)＝P ^-(k)M(k) ^T(M(k)P ^-(k)M(k) ^T+R) ^-1 (16)

4. correcting state is estimated

$\hat{x}\left(k\right)={\hat{x}}^{-}\left(k\right)+K\left(k\right)(y\left(k\right)-{C^{\&prime;}}_d{\hat{x}}^{-}\left(k\right))---\left(17\right)$

5. upgrade error covariance matrix P (k):

P(k)＝(I-K(k)M(k))P ^-(k) (18)

Wherein

$G\left(k\right)={\frac{\&PartialD;f}{\&PartialD;x}|}_{x=\hat{x}(k-1)}$

$=\left[\begin{array}{ccc}1-\frac{T_s}{{\mathrm{\&tau;}}_r}& -{\widehat{\mathrm{\&omega;}}}_r(k-1)T_s& -{\widehat{\mathrm{\&Psi;}}}_{\mathrm{\&beta;r}}(k-1)T_s\\ {\widehat{\mathrm{\&omega;}}}_r(k-1)T_s& 1-\frac{T_s}{{\mathrm{\&tau;}}_r}& {\widehat{\mathrm{\&Psi;}}}_{\mathrm{\&alpha;r}}(k-1)T_s\\ 0& 0& 1\end{array}\right]---\left(19\right)$

$M\left(k\right)={\frac{\&PartialD;h}{\&PartialD;x}|}_{x={\hat{x}}^{-}\left(k\right)}=\left[\begin{array}{ccc}-\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& -{{\widehat{\mathrm{\&omega;}}}^{-}}_r\left(k\right)\frac{L_m}{L_r}& -{{\widehat{\mathrm{\&Psi;}}}^{-}}_{\mathrm{\&beta;r}}\left(k\right)\frac{L_m}{L_r}\\ {{\widehat{\mathrm{\&omega;}}}^{-}}_r\left(k\right)\frac{L_m}{L_r}& -\frac{L_m}{L_r{\mathrm{\&tau;}}_r}& {{\widehat{\mathrm{\&Psi;}}}^{-}}_{\mathrm{\&alpha;r}}\left(k\right)\frac{L_m}{L_r}\end{array}\right]---\left(20\right)$

Subscript "-" among formula (13)-Shi (20) is all represented premeasuring, and subscript " ∧ " is all represented estimator; Variable, meaning of parameters are all identical with formula (1);

So far, obtained asynchronous machine depression of order EKF rotating speed algorithm for estimating at each sampling period T _sInterior concrete calculation procedure and computing formula are for the hardware language description of next step FPGA is got ready;

In the 3rd step, design FPGA realizes the algorithm structure that depression of order EKF rotating speed is estimated:

According to formula (14), control u (k) multiply by B ' _dBattle array (96) back and state transitions battle array A ' _d(k) and state

Long-pending (97) obtain the result that 1. EKF the goes on foot status predication through adder (95) addition

Wherein, u (k) is the stator current vector; Wherein

It is the state estimation that a last bat obtains

Obtain through unit delay unit (98) delay one bat; Predicted state

With a measurement battle array C ' _d(k) multiply each other (91), subtracter (92) be used for calculating measuring amount y (k) therewith the difference of product give Kalman gain adjustment module (93), i.e. the of EKF 3. step; Before this, FPGA should trigger the 1. the step calculate in the computing circuit of trigger error Estimates on Covariance Matrix, upgrade P by formula (15) ^-(k), promptly accomplish the 2. step of EKF; Use Kalman gain adjustment module (93) to predicted state

Proofread and correct by formula (17), adder (94) is used on predicted state, adding the correcting value by the long-pending decision of Kalman gain and measure error, obtains final state estimation

Promptly accomplish the 4. step of EKF; After this, the computing circuit that FPGA answers the trigger error covariance matrix to upgrade upgrades P (k) by formula (18), promptly accomplishes the 5. step of EKF; So far, all interior calculating of sampling period of EKF have been accomplished;

The computing formula of measuring amount y (k) is following:

$y\left(k\right)=\left[\begin{array}{c}{u}_{\mathrm{\&alpha;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&alpha;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&alpha;s}}\left(k\right)\\ u_{\mathrm{\&beta;s}}\left(k\right)-(R_s+\frac{{L_m}^2}{L_r{\mathrm{\&tau;}}_r})i_{\mathrm{\&beta;s}}\left(k\right)-\mathrm{\&sigma;}{L}_s{\stackrel{\&CenterDot;}{i}}_{\mathrm{\&beta;s}}\left(k\right)\end{array}\right]---\left(21\right)$

Plus and minus calculation in subtracter (991) and subtracter (992) perfect (21), stator current measured value i _S(k)=[i _{α s}(k), i _{β s}(k)] ^TMultiply each other obtain long-pending and the differential signal of the stator current that calculates by weighted network (994), coefficient (996) with coefficient (993) With the product of coefficient (995) as subtrahend and stator voltage command value

Poor, measuring amount y (k) just; Differential signal to stator current Calculating adopt method---the three point value differential of evaluation differential in the numerical analysis:

${\stackrel{\&CenterDot;}{i}}_s\left(k\right)=\frac{1}{2T_s}(i_s(k-2)-4i_s(k-1)+3i_s\left(k\right))---\left(22\right)$

Displacement weighted sum in the following formula (22) is accomplished i by weighted network (994) _S(k) each is clapped and is postponed to realize that by a unit delay unit this weighted sum and coefficient (996) multiply each other and be the differential of stator current

The 4th step,

Based on the algorithm structure of above description, FPGA is carried out hardware language VHDL describes, the equal using system clock of all modules CLK is as clock signal, Rstx signal that each module receives and Cwel signal deciding their calculating sequential; The Rstx signal is low effective enable signal, and Cwel is a register type, count range 0 to 12, and the duration of each count value is 3 CLK, the matrix multiplication that is designed can be accomplished at 3 CLK; In each count cycle; When Rst1 when low, measuring amount computing module (937), status predication module (932) and error covariance matrix estimation module (936) are enabled, computing begins; The operation time that error covariance matrix estimation module (936) needs is the longest, is 18 CLK; After 18 CLK; Be that the Cwel count value is after 6; 1. the of the calculating measuring amount computing module (937) of measuring amount y (k) and EKF go on foot status predication module (932) and the and 2. go on foot error covariance matrix estimation module (936) and all carried out computing concomitantly, for the of EKF 3. go on foot Kalman gain calculating (935) all set when time each variate-value; The Cwel count value is 7～9 these 3 count values, and Rst3 is when low, and Kalman gain calculating (935) is enabled; The Cwel count value is 10～12 these 3 count values; And when Rst4 is low; State correction module (933) and error covariance matrix update module (934) are enabled simultaneously; The Cwel count value is after 12, and 4. of EKF goes on foot state correction module (933) and and 5. go on foot error covariance matrix update module (934) and all carried out computing concomitantly, as inferior state estimation value i _{α s}, i _{β s}, Ψ _{α r}, Ψ _{β r}, ω _rBe ready to be sent to the parallel port and passed master control DSP back; i _{α s}, i _{β s}Be the estimated value of α, β axle stator current, Ψ _{α r}, Ψ _{β r}Be the estimated value of α, β axle rotor flux, ω _rEstimated value for rotor electrical angle speed.

Images