MEASUREMENT OF TRANSFER FUNCTION IN A MECHATRONIC SYSTEM
20230082667 · 2023-03-16
Inventors
Cpc classification
G05B19/4155
PHYSICS
International classification
Abstract
A device comprising at least one electromechanical system and an electronic control unit comprising both a servo control loop connecting together the input and the output of the electromechanical system, and also a digital corrector. The device includes a software module having an input and an output connected respectively to the output and to the input of the electromechanical system. For each frequency of a set of frequencies that are to be analyzed, the software module is programmed: to issue a sinusoidal excitation signal on its output and to recover a response signal on its input; to model the response signal by means of a Fresnel representation and to filter it about the excitation frequency by means of at least one Kalman filter in order to obtain a state vector of the filtered response signal; and to obtain at least one transfer function for the system from the model, from the excitation signal, and from the input signal.
Claims
1. A device comprising at least one electromechanical system having an input and at least one output, and an electronic control unit comprising both a servo control loop connecting together the input and the output of the electromechanical system, and also a digital corrector, the device being characterized in that it includes a software module having an input connected to the output of the electromechanical system and an output connected to an input of electromechanical system, and in that, for each frequency in a set of frequencies that are to be analyzed, the software module is programmed: to generate at least one sinusoidal excitation signal (P); to issue the sinusoidal excitation signal on its output and to recover a response signal on its input; to model the response signal by means of a Fresnel representation and to filter it about the excitation frequency by means of at least one Kalman filter in order to obtain a state vector of the filtered response signal; to multiply the state vector by a rotation matrix in order to obtain a real part and an imaginary part; and to perform complex number division to obtain at least one transfer function for the system from the model, from the excitation signal, and from the input signal.
2. The device according to claim 1, wherein the signal generator is based on a representation of a Fresnel vector rotating in compliance with the following model: initial condition:
X.sub.k+1.sup.P=ΔA.sub.pX.sub.k.sup.P the signal to be applied to the output of the corrector:
P.sub.k=C.sub.P.Math.X.sub.k.sup.P with the following matrices for the state representation:
3. The device according to claim 1, wherein the dynamic range of the excitation signals and the dynamic range of the measurement signals are defined by the following equations:
4. The device according to claim 3, wherein the matrices ΔA.sub.M, B.sub.M, C.sub.M, D.sub.M are determined as follows:
5. The device according to claim 3, wherein the matrices ΔA.sub.M, B.sub.M, C.sub.M, D.sub.M are determined as follows:
6. The device according to claim 3, wherein the matrices ΔA.sub.M, B.sub.M, C.sub.M, D.sub.M are determined as follows:
7. The device according to claim 3, wherein the matrices ΔA.sub.M, B.sub.M, C.sub.M, D.sub.M are determined as follows:
8. The device according to claim 3, wherein the Kalman filter is implemented in the following simplified form:
9. The device according to claim 2, wherein the covariance matrix of the measurement noise is selected to be equal to 1 and the covariance matrix W of the model noise is selected as follows:
10. The device according to claim 9, wherein the variance σ.sub.bias.sup.2 is set at 10.sup.−2 and the variance σ.sub.line.sup.2 is set at 10.sup.−4 for a sampling frequency equal to 3 kHz.
11. The device according to claim 2, wherein the rotation matrix is defined as a function of the components of the vector that is representative of the sinusoidal excitation signal.
12. The device according to claim 11, wherein the rotation matrix has the following form:
Description
BRIEF DESCRIPTION OF THE DRAWINGS
[0024] Reference is made to the accompanying drawings, in which:
[0025]
[0026]
[0027]
[0028]
[0029]
DETAILED DESCRIPTION OF THE INVENTION
[0030] With reference in particular to
[0031] The motors 3 and 4, and the associated angle sensors 3 bis and 4 bis are connected to a control device 10 that issues respective commands m1 and m2 to the motors 3 and 4, and that receives respective measurements c1 and c2 from the angle sensors 3 bis and 4 bis.
[0032] The optical sensor 1 is provided with at least one inertial angle sensor 5 in order to detect an angular velocity c3 of the aiming axis v relative to a first reference axis that is collinear with the elevation axis y, and an angular velocity c4 of the aiming axis v relative to a second reference axis that is perpendicular to the aiming axis v and to the elevation axis y. In this example, the inertial angle sensor is a gyro having two sensing axes, such as the gyro produced under the trademark GSL by the Applicant. In a variant, any other angle measuring device could be used, and for example it is possible to use two inertial sensors, each having a respective sensing axis, such as the device produced under the trademark QUAPASON by the Applicant. The inertial sensor could be a mechanical sensor having a vibrating resonator, e.g. a hemispherical resonator, or use could be made of other technologies, such as optical fibers (a so-called fiber-optic gyro (FOG) sensor).
[0033] The inertial sensor 5 is also connected to the control device 10, and it transmits the measurement signals c3 and c4 thereto.
[0034] With reference also to
[0035] For each linear system referenced H, the control device 10 includes an electronic control unit 11 having a servo control loop that connects together the input and the output of the linear system by means of a digital corrector 12. The electronic control unit 11 is thus arranged to act via the digital corrector 12 to control the motors as a function of commands coming from a user interface 20 (visible in
m=K*(setpoint−c) and c=H*m.
[0036] This mode of operation is conventional and is not described in greater detail herein.
[0037] The control device 10 exchanges data with memory modules 30 and 40. The memory module 30 comprises a table of the frequencies f.sub.p that are to be scanned, of the amplitudes P.sub.0 of the associated excitation signals, and of the vectors of the associated Kalman gains K.sub.f.
[0038] The control device 10 includes a software module 13 having an output p connected to an output e3 of the digital corrector 12, at least one first measurement input e1 connected upstream from an input of the digital corrector 12, and two measurement inputs e2 and e3 connected to the output of the digital corrector 12 respectively downstream and upstream from the connection between said output of the digital corrector 12 and the output p of the software module 13. Thus, the output p added to the output e3 of the digital corrector 12 gives the command e2 that is to be applied to the motor via the digital-to-analog converter and the signal m.
[0039] The software module 13 is programmed: [0040] to form an excitation generator based on a rotation increments matrix and on a Fresnel model, the generator being arranged to issue (to an input p of the system S) final excitation signals P of different frequencies and/or of different amplitudes on the output of the corrector by using a software sum e2=e3+P; [0041] to recover on each of its inputs a respective sinusoidal measurement signal e1 to e6 for each of the sinusoidal excitation signals; [0042] to model the measurement signals e1 to e6 (
[0046] The software module 13 is arranged to scan a plurality of frequencies, recovering for each frequency the excitation amplitude and the gains of the Kalman filters stored in a memory (30 see
[0047] By way of example, the software module 13 is arranged to be capable of estimating a plurality of transfer functions, which are identified as FT1 to FT6 in
[0048] The principle of the invention for determining one of the transfer functions is described below in its three main stages: generating a sinusoidal excitation signal at the frequency f.sub.p, bandpass filtering using a Kalman filter about the frequency f.sub.p, and multiplying by a rotation matrix in order to obtain the real and imaginary parts of the transfer function between the excitation and the measurement under consideration at the frequency f.sub.p.
[0049] The excitation signal (disturbance) is a sinusoidal signal, for example a cosine of amplitude P.sub.0 at the frequency f.sub.p, and it can be modelled in the form of a vector that is representative of a complex number and that rotates at a constant speed {dot over (θ)}=2π.Math.f.sub.P (known as a “Fresnel representation”). At the instant t.sub.k, the complex number is written as follows:
{circumflex over (P)}.sub.kP.sub.0.Math.exp(i.Math.2π.Math.f.sub.P.Math.t.sub.k)=P.sub.0.Math.exp(i.Math.θ.sub.k)=P.sub.0.Math.cos(θ.sub.k)+i.Math.P.sub.0.Math.sin(θ.sub.k)
[0050] The voltage disturbance that is actually injected may be expressed as being the real part of the complex number, i.e.:
P.sub.k=P(t.sub.k)=P.sub.0.Math.cos(θ.sub.k)=Re{P.sub.0.Math.exp(i.Math.2π.Math.f.sub.P.Math.t.sub.k)}
[0051] While the imaginary part of the complex number is left free:
P.sub.0.Math.sin(θ.sub.k)=Im{P.sub.0.Math.exp(i.Math.2π.Math.f.sub.P.Math.t.sub.k)}
[0052] This can be written in matrix form as follows:
[0053] For a constant frequency f.sub.p, the transition matrix between two sampling periods T.sub.e corresponds to rotation through the angle Δθ=θ.sub.k+1−θ.sub.k=2π.Math.f.sub.P.Math.T.sub.e, where T.sub.e is the sampling period of the calculation means (T.sub.e=t.sub.k+1−t.sub.k). The transition matrix is written as follows:
[0054] The matrix relationship between two consecutive instants is thus:
[0055] Below, it should be understood that the rotation matrix at the instant k is the product of k transition increments:
[0056] It should be observed that the inverse rotation matrix is simply the transpose. It too is calculated using the components of the excitation vector:
[0057] The notation (.circle-solid.).sub.n designates the “n.sup.th” component of the vector representative of the disturbance.
[0058] Taking the disturbance P.sub.k as the first component of the state vector, the following discrete dynamic system is obtained:
in which B=D=0 and C=[1 0].
[0059] As shown in
[0060] The matrix form of the dynamic system is adapted for performing Kalman filtering. P.sub.0 is taken as the initial amplitude condition of the excitation signal (disturbance) associated with a phase of zero, such that the initial vector is:
[0061] Concerning the frequency response of the electromechanical system and the matrix product, a formulation in the form of a linear dynamic system makes it possible to treat the physical system as a linear filter that receives the disturbance P.sub.k as input to give the measurement M.sub.k as output. Using the z transform and the notation {circumflex over (F)}(z) for the transfer function that is to be estimated, the following is obtained:
{circumflex over (M)}(z)={circumflex over (F)}(z){circumflex over (P)}(z)
[0062] When the delay operator z.sup.−1=exp(−i.Math.2π.Math.f.Math.T.sub.e)∀f, this gives:
M.sub.k−1=z.sup.−1.Math.M.sub.k
[0063] For a linear system excited with a pure cosine at the frequency f.sub.p, the measurement is a cosine of amplitude M.sub.P and of phase ψ.sub.P, and it takes the following matrix form:
[0064] The measurement that is directly accessible is M.sub.k=C.Math.X.sub.k.sup.M, whereas it is the following that is being looked for:
[0065] This can be written in matrix form:
[0066] Where the “matrix” measurement X.sub.k.sup.M and the transfer function X.sup.F are associated by the following trigonometrical relationships (see the Fresnel diagram of
[0067] I.e., on being inversed:
[0068] The real and imaginary components of the transfer function are obtained merely by a matrix product of the ideal measurement signal multiplied by a single rotation matrix of components that have already been calculated by generating the excitation signal:
{circumflex over (F)}(f.sub.P)=(X.sup.F).sub.1+i(X.sup.F).sub.2
where i is the imaginary unit such that i.sup.2=−1.
[0069] There follows a description of the bandpass filtering by Kalman filtering.
[0070] In order to estimate the real part and the imaginary part of the measurement X.sub.k.sup.M from the real measurement M.sub.k, use is made of a converged-gain Kalman filter that acts as a selective bandpass filter about the excitation frequency:
[0071] In which: [0072] w.sub.k is the noise of the scalar model; [0073] v.sub.k is the scalar measurement noise; [0074] ΔA is the transition matrix; [0075] C=[1 0] is the observation matrix of the disturbance generator; and [0076] B.sup.T=[1 1] is the state noise injection matrix.
[0077] The variance of the measurement noise is V=E[v.sub.k.sup.2].
[0078] The variance matrix of the model noise is:
W=BE[w.sub.k.sup.2]B.sup.T
[0079] The filter may be represented in the form of the following state representation:
[0080] The matrix K.sub.f is the converged-gain matrix of the Kalman filter, such that:
K.sub.f=ZC.sup.T(CZC.sup.T+V).sup.−1
and
ΔAZ(ΔA).sup.T−Z−ΔAZC.sup.T(CZC.sup.T+V).sup.−1CZ(ΔA).sup.T+W=0
[0081] In which Z is the solution of the Riccati equation.
[0082] In practice, in order to determine the transfer function of the linear system in real time, the following are known, the frequency f.sub.P of each disturbance signal, the amplitude P.sub.0, the duration of the filtering T.sub.f, and the gain of the Kalman filter, i.e.:
[0083] In order to scan all of the frequencies in a given interval, the following algorithm is performed:
TABLE-US-00001 Algorithm: Review all of the frequencies f.sub.P • If t.sub.k = 0
Output
[0084] X.sup.F=[real FT,imag FT].sup.T for each frequency and for each measurement.
[0085] It should be observed that in this example, the equations of the Kalman filter are put in the form of a recursive filter with a single state vector, as in equation 2.31 of the document “Introduction to Kalman filtering” by D. Alazard. Nevertheless, it is also possible to write them in the prediction/updating form with two state vectors (see equations 2.23 and 2.25 of that document) as in the description below of the continuous component.
[0086] Estimates are thus obtained of at least three transfer functions of the looped system, between the excitation signal and at least three measurements, such that
without any need for operator intervention in order to separate the servo control loop from the actuator.
[0087] The various transfer functions of the loop are then determined as follows.
[0088] With the notation of
a/ For the current frequency, storing for the current frequency the real parts of at least the following three closed loops:
b/ Calculating the open loop with:
c/ Calculating the transfer function of the electromechanical portion:
[0089] There follows a description of a preferred embodiment, making use of reformulation and extension in order to take the continuous components of the measurements into account.
[0090] In the preferred embodiment, it is not a cosine that is injected, but rather an excitation sine having a predetermined frequency f.sub.P in hertz (Hz).
[0091] When a sine type excitation signal is injected into a linear system, the system responds with a sine line that is amplified by a gain
and phase shifted by a phase Ψ.sub.p.
[0092] The state representation of the sinusoidal line of the measurement without a continuous component thus has two dimensions and is given by:
where
[0093] The state vector represents the cosine (plotted along the x-axis) and the sine (plotted along the y-axis) on a circle of radius M.sub.p=G.sub.pP.sub.0 for a vector rotating at the frequency f.sub.p and starting from an initial position determined by the angle Ψ.sub.p.
[0094] Thus, the following matrices are obtained for the state representation of the Fresnel vector:
where, in conventional manner for state representation, ΔA is the transition matrix, B is the control matrix, C is the observation matrix, and D is the direct action matrix.
[0095] Nevertheless, in practice, it is probable that the response of the system also includes a continuous component corresponding to the operating point about which the system is excited.
[0096] This continuous component constitutes a bias that has the following discrete representation:
[0097] By adding this representation to the state representation of the sine line, the state representation is obtained for a sine having a continuous component in three-dimension for the various measurements, namely:
where:
[0098] There is no need for the excitation generator to have a continuous component. In order to design the sinusoidal excitation generator, it suffices to simulate the model described by the state representation of the sine line. The amplitude and the phase are selected by imposing initial conditions on the model. A phase of zero is selected for the sine wave that is to be generated so as to start excitation of the system progressively (sin(t=0)=0). A generator of a sine wave of amplitude P.sub.0 is obtained as follows: [0099] initial condition
X.sub.k+1.sup.P=ΔA.sub.PX.sub.k.sup.P
P.sub.k=C.sub.PX.sub.k.sup.P [0101] and with the following state representation matrices:
[0102] Each measurement from a sensor is a sine line PGP-including a continuous component, from which the state of the Fresnel vector is estimated by using a Kalman filter constituting an estimator corresponding to a bandpass filter about the frequency f.sub.p for the first two components of the state vector, and a lowpass filter for the continuous portion, based on the models (1) and (2).
[0103] Once this state has been estimated, it is possible to estimate the corresponding transfer function between the excitation (input) and the measurement (output).
[0104] At the output from the transfer function estimator, there are obtained directly in the real part and the imaginary part of the transfer function at the frequency f.sub.p.
[0105] The model of the augmented Kalman estimator takes the following general form:
[0106] With: [0107] ΔA.sub.M, B.sub.M, C.sub.M, D.sub.M as above; [0108] U.sub.k represents the commands; [0109] v.sub.k the measurement noise; [0110] w.sub.k the noise of the model; and [0111] the selected matrix M is the following:
[0112] Since the matrices B.sub.M, D.sub.M are zero, it is possible to simplify the model:
[0113] The notation V[1×1] and W[3×3] is used for the covariance matrices of the measurement noise v.sub.k and of the model noise w.sub.k, which noise is assumed to be Gaussian for a linear Kalman estimator. The steady state of the estimator is determined solely by the selected covariance matrices V and W. Given that the estimator is implemented with converged gain, it is possible to set V=1 arbitrarily without loss of generality. W, the covariance matrix of the noise of the model, is selected as follows:
[0114] Thus, selecting the variance pair (σ.sub.line.sup.2 σ.sub.bias.sup.2) serves to define completely the steady state of the Kalman estimator.
[0115] A state representation of the estimator is written:
[0116] The following is then obtained:
where:
K.sub.f=ZC.sub.M.sup.T(C.sub.MZC.sub.M.sup.T+V).sup.−1
[0117] In this representation, Z is the only solution that stabilizes the algebraic Riccati equation:
Z=ΔAZ(ΔA).sup.T−ΔAZC.sub.M.sup.T(C.sub.MZC.sub.M.sup.T+V).sup.−1C.sub.MZ(ΔA).sup.T+W
[0118] It should be observed that this equation presents a singularity in the event that the frequency f.sub.p is equal to half the sampling frequency F.sub.e=1/T.sub.e. In the general case, this equation can be solved for example using the “dare” function of the “Matlab” language that also serves to determine the eigenvalues of A.sub.k that enable the time constant τ.sub.K of the bandpass filter to be determined for the first two components of the state vector and of the lowpass filter of the continuous component.
[0119] As mentioned above, the settings of the estimator (i.e. the width of the band about f.sub.p, and the cut-off frequency for the lowpass filter) are obtained by selecting the pair of variances (σ.sub.line.sup.2 σ.sub.bias.sup.2) that enable the steady state Kalman estimator to be defined completely.
[0120] Assuming that the variance σ.sub.bias.sup.2 is fixed, then reducing the variance σ.sub.line.sup.2 gives rise to a reduction in the width of the bandpass filtering, and thus to a better estimate of the sine line as a result of better filtering of the noise and of any nonlinearities of the system. If the time constant τ.sub.K of the filter is equal to the time constant of the line, then this reduction in the variance σ.sub.line.sup.2 is accompanied by an increase in the time constant of the filter. It is therefore necessary to seek a compromise between the quality of the estimate (narrowness of the filter) and the convergence time of the estimator.
[0121] Assuming that the variance σ.sub.line.sup.2 is fixed, then a reduction in the variance σ.sub.bias.sup.2 give rise to a deterioration in the filtering of low frequencies for the cosine and sine components. In addition, the cut-off frequency of the response corresponding to the continuous component (bias) decreases when the variance σ.sub.bias.sup.2 decreases (which amounts to increasing measurement noise relative to the noise of the model), giving rise to filtering that is more conservative.
[0122] By way of example, for a given sampling frequency F.sub.e, for given minimum and maximum frequencies f.sub.min and f.sub.max of the excitation signal, and for a number N.sub.points of frequencies to be identified in this frequency range, it is proposed to set the estimator as follows: [0123] determine the variance σ.sub.line.sup.2 so as to obtain the desired width for the bandpass filtering. This is done more for middle/high frequencies in the frequency range so that the filter is not disturbed by the continuous component of the model. It is ensured that the filter is indeed symmetrical about the selected frequency; [0124] determine the variance σ.sub.bias.sup.2 on the lowest frequency that is to be identified (f.sub.min) so that the selectivity of the filter remains acceptable at this frequency; and [0125] calculate the total duration for the identification induced by these settings. If the generation is not satisfactory, increase the variance σ.sub.line.sup.2 or decrease the number of frequencies to be identified. This type of scanning has a strong impact on the duration of identification: in order to reduce it, it is preferable to select scanning that is linear rather than logarithmic.
[0126] In the present example, (σ.sub.line.sup.2 σ.sub.bias.sup.2)=(10.sup.−4 10.sup.−2) is selected for a sampling frequency F.sub.e of 3 kHz. By way of example, this setting is appropriate for aiming an optronic device. For some other application, the above-mentioned values can serve as starting points for converging on settings that are more appropriate for said application.
[0127] In particular, such an application might require a sampling frequency that is different. However, it is known that for a given setting of the pair (σ.sub.line.sup.2 σ.sub.bias.sup.2), the width of the bandpass filter increases with increasing sampling frequency: identification is thus faster but of poorer quality. This problem can be solved by recalculating the values of the (σ.sub.line.sup.2 σ.sub.bias.sup.2) pairs so as to guarantee a time constant that is almost identical and a frequency response that is equivalent when changing from one sampling frequency to another.
[0128] In a variant, it is possible to reduce the calculation burden by making use of an estimator that is no longer in the form of a state representation
as described above (i.e. in the form of a recursive filter having a single state vector, but rather in the prediction/updating form with two state vectors, by using the notation {circumflex over (X)}.sub.k+1|k.sup.M for the predicted state and {circumflex over (X)}.sub.k+1|k+1.sup.M for the updated state (output from the estimator), with the recurrence equations of the filter at instant k+1 being:
[0129] The Kalman gain K.sub.f is the same as that defined above. This implementation is mathematically equivalent to the preceding implementation, but less expensive in calculation time.
[0130] The transfer function estimator calculates the real and imaginary parts of the transfer function {circumflex over (F)}, between the sinusoidal excitation signal and the measurement that is obtained, for one of the following three closed loops:
[0131] It is thus sought to estimate the following (where i is the imaginary unit such that i.sup.2=−1):
{circumflex over (F)}(f.sub.p)=R.sub.p+i I.sub.p=G.sub.p(cos(Ψ.sub.p)+i sin(Ψ.sub.p))
[0132] The state of the excitation generator at instant k is written X.sub.k.sup.P, and the process estimated by the Kalman filter on the basis of the measurement is written {circumflex over (X)}.sub.k|k.sup.M. The following is obtained:
[0133] The pair
is obtained from the state of the process {circumflex over (X)}.sub.k|k.sup.M by a matrix T constructed using the terms X.sub.k.sup.P and P.sub.0:
[0134] As before;
is a rotation matrix.
[0135] This gives:
giving rise to estimates of the real part R.sub.p and of the imaginary part I.sub.p such that:
[0136] This estimation should not be performed until the Kalman estimator has converged sufficiently. In order to be sure that it has, a length of time is allowed to elapse that is equal to several times the characteristic time constant τ.sub.K of the estimator before estimating the transfer function. By taking a length of time that is equal to six times the characteristic time constant τ.sub.K, it is ensured that the estimation error has convergence of 0.25%.
[0137] In order to ensure that the bias of the mechanical system (e.g. associated with the presence of unbalance or of a return member in a movement transmission) does not disturb determination the transfer function, it should be observed that it is necessary to ensure that the amplitude of the sinusoidal signal is sufficient to be able to ignore the effects that are not generated by the sinusoidal signal.
[0138] With reference to
[0139] Naturally, the invention is not limited to the implementation described, but covers any variant coming within the ambit of the invention as defined by the claims.
[0140] In particular, the invention is applicable to any type of electromechanical system that can be modelled by means of a linear transfer function, and it is therefore not limited to an application to an airplane control surface or to a device for aiming and stabilizing an optical sensor. In particular, the invention can be used with any device including an actuator controlled by servocontrol and measurements that are accessible to calculation means.
[0141] Other modelling parameters may be selected as a function of the intended application. The sampling frequency, the variances, the number of frequencies to be estimated, . . . .
[0142] The term “sinusoidal signal” is used to cover both a sine signal and a cosine signal.
[0143] The software module may receive as input the measurement of the output variable of the linear system and/or the difference between the measurement of the output variable of the linear system and the control setpoint.
[0144] Although in the embodiments described the rotation matrix is defined as a function of the components of the vector representative of the sinusoidal excitation signal, which is advantageous in terms of simplifying calculation, that is not essential.
[0145] Furthermore, it should be observed that in the embodiment described, account is taken of three measurements e1, e2, and e3: [0146] e2 and e3 for the “open loop” (K*H), i.e. one upstream from the summing circuit (e3) and one downstream from the summing circuit (e2); and [0147] e1 and e2 for the transfer function of the electromechanical system (H).
[0148] In contrast, if it is desired to measure simultaneously both inertial speeds, the positions of the two coders, and both motor currents, it is possible to use e4, e5, and e6 in addition to e1.
[0149] The matrix T may have various different forms: [0150] to have b.sub.0, use: