Team,
Presently I am doing work on multicopter modeling and control. One topic of interest to me is how can the tilt dynamic model of a multicopter be accurately determined from flight data, without disabling the controls or injecting any sort of intrusive signals. It turns out it can be done from normal flight data using a few basic concepts from signal processing. The method is described here. The above plot is the measured pitch dynamic impulse response of my draganflier. The x axis is time. The y axis is the pitch rate impulse response to an impulse command sent to the ESC. There are two curves in the plot. h(t) is the response, computed from flight data. H(t) is an exponentially decaying sinusoid fitted to the measurement. In this case, the Laplace transform of the impulse response is 5/(s*s+4*s+29), which is of exactly the form you would expect to arise from an analysis of fixed pitch control. Basically, in the case of my draganflier, there is about a 0.25 second delay between the time a new tilt rate is commanded to the motors and when quad begins tilting. I have measured a similar delay for the ArduCopter. The cause of this delay is mainly the inertia of the propellers, which delays the response of propeller thrust to change in motor torque, and the inertia of the quad, which delays the tilt rate response of the quad to propeller thrust.
As a check on the accuracy of the method, it is possible to convolve H(t) with control inputs to predict tilt rate. The following is a comparison of measured and predicted pitch rate during a flight of my draganflier.
Don't get too excited about the way the model seems to predict the high frequency features. If you are curious about that, read my report.
From the form of the tilt dynamic model it is possible to determine the form of an optimal feedback controller. It is not a PID controller. But that will be the subject of another report.
Best regards,
Bill
Comments
The different vector matrix methods are a neat discussion. What I have found important is that there needs to be little mechanical vibration to the IMU.
Hi William and hi all,
I think you can get even more accurate results using a similar strategy but working in the Frequency Domain. Here is what you can do (without so much effort):
1. Conduct experiments where you give inputs that correspond to a frequency sweep. You can add this to your input signal (add a chirp signal there) but even if you do not want to do this then at least you should be able to provide the needed frequencies through your reference commands. The controller is going to mix them through the feedback of the data but more or less they are going to be there. As long as the order of the model is higher than 1 and you are recording the control actions then it is possible to identify the model.
2. Try to conduct the experiments in a way that you start and end at trim.
3. Use frequency-domain based tools to check if your excitation was good enough.
4. Develop a Grey-box of your model. This means that you take the 6DOF model or decoupled verions of it (Roll - Pitch - Yaw and so on). So you have a specified structure but unknown parameters.
5. Use a regressor (MATLAB is good for that - it also has a special tool for system id) and solve the problem in the FREQUENCY DOMAIN. This means that you solve the regression between the frequency response of the experimental flight response and the estimated frequency response of the model to be identified. Frequency domain poses several advantages, and the first (and very related with the fact that you are working with feedback enabled) is the frequency seperation and the capability to focues on the frequencies produced from the vehicle and not from noise.
6. Validate the model in time-domain as well as frequency domain.
The aforementioned general steps are ok even for the open-loop unstable model, which is another key feature of the frequency domain.
Regarding the online implementation you have to look for an ARX model. Keep the model order low and you can generally ensure convergence of the parameters. The open-loop unstable dynamics (critically stable) are going to be quite tricky then. Are you also interested for adaptive control tuning?
I am not trying to attenuate anything in the signal
But effectively you are. And theoretically you have to in order to prevent aliasing in you 1/40 decimation stage.
For what it is worth, I suppose the modern heli tail gyros as well as the flybarless 3-axis gyro systems available today, are using this hi frequency sampling together with integration (just like William has it done). This might explain why those gyro systems manage to maintain very low gyro drift while installed in harsh vibration environments (nitro/ gas) engines on copters performing violent 3D-flight manouvers (3D-flying introduces vibrations due to brutal cyclic usage of the rotor disc, and flybarless rotor heads are quite rigid). Without help of secondary reference like magnetometer or accelerometers.
To me this way to deal with sampling and integration seems very sound. I think the sooner it could be introduced into ArduCopter realm, the better.
Cheers / Tomas
Bill,
we also measure drift in degrees/s, what I am talking about is the slope of omegaI. The y axis of the graph I showed you is in degrees/s, but its average slope is around 0.02 degrees/s/s. My apologies for not labelling my axes! As I understand it, the ki value in DCM that controls the growth of omegaI is essentially a limit on that slope.
Regarding the effect of the mag on omegaIz, I agree completely, and it is a problem I'm trying to deal with right now in the context of the Madgwick quaternion system (I have a compile time option to switch between DCM and the Madgwick MARG system).
Quite possibly. I rewrote our gyro calibration code yesterday because I suspected that it was producing a poor value. I'm going to watch the logs from users to see if the initial slope of omegaI decreases as a result. I think it will.
In APM2 we use a Invensense MPU6k. In APM1 we use a IDG500 for x/y and a ISZ500 for z.
I don't think it is temperature effects however. In our simulator I've been able to produce very similar 'drift' just by injecting white noise. I suspect the PI values we're using to control omegaI are not good and our 1kHz sample rate is not high enough to control it, combined with a base drift from poor calibration.
We've only just started recording the chip voltage during flight, so we don't have much data yet. I'm hoping that future user logs will allow us to see if the logged chip voltage is in any way correlated with the accumulation of error in the DCM (we're now logging the magnitude of the errors terms in DCM, as well as the omegaI integrator components).
Cheers, Tridge
Hi Kari,
Yes, I am using a 200 tap mean filter. Just to be clear, that is the theoretically correct thing to do for the critical downstream computations that are basically integrators. I am not trying to attenuate anything in the signal, I am treating everything as signal, and then integrate it. That has provided incredible performance.
And yes, the data in the spreadsheet is at 40 Hz.
Best regards,
Bill
O.k. so basically you are using 200 tap mean filter for the decimation.
That filter gives you approx -45dB attenuation in the stop band. You can actually do much better with 200 taps filter, easily reaching -80dB attenuation in the stop band.
The raw data in the spreadsheet (pfb column) is it at the raw 8000Hz sampling rate or is it after the decimation stage?
Andrew,
One more observation. Part of the problem could be the reference vectors themselves. Non-zero values of omegaI simply mean that the gyros and the reference vectors disagree.
For example, I find that uncompensated yaw drift is essentially zero. But if I hook up a magnetometer, I will see a non-zero value of omegaIz, because of magnetometer errors.
So, the problem that we face is that the gyros and reference vectors have different sorts of errors, and once you rely on a reference vector for drift compensation, you are stuck with the errors it introduces.
In your case, I suspect that accelerometer errors are contributing to what you see in roll and pitch, and magnetometer errors are contributing to yaw.
When I have some time, I could look into this on my end with UDB and MatrixPilot, so that we could compare notes. Right at the moment, though, I am tied up with a bunch of other things, including some things that I am working on for GE. (They have rehired me back out of retirement.)
Best regards,
Bill
Hi Andrew,
You and I were responding at the same time.
Just to be clear, I am running without any yaw drift compensation. So, my omegaIZ is identically zero. I am running without a GPS and without a magnetometer. Even so, MatrixPilotQuad will hold yaw orientation to within a couple of degrees for 10 minutes, and land in the same orientation as takeoff.
One thing you said really caught my attention...you are measuring gyro drift rates in units of degrees/sec/sec, where I measure them in units of degrees/sec. Which suggests that the drift mechanisms that you are seeing and the ones that I am seeing are much different.
In my case I think that it is mostly that the recording of the initial offsets are ever so slightly off, because of noise.
In your case I am not sure what is going on. Temperature effects maybe? What kind of gyros are you using?
Finally, I have not looked at the values of omegaI in MatrixPilot lately, I will put that on my list of things to do. Last time I looked, they were rather small.
By the way, one of the things we found to be the biggest source of gyro drift with the UDB is power supply voltage. We use a 3.3 volt regulator for the gyros and accelerometer. However, some setups introduce a large amount of noise on the 5 volt bus that spills down onto the 3.3 volt bus. This causes major interference with gyro drift. So, the UDB4 offers an option to split the 5 volt bus into two sections, one for connections to the radio, and another to the connections to the servos. Since you are not using servos on a quad, I do not suspect that is a factor in your case. But its something to think about.
Best regards,
Bill
Team,
Thank you all so much for your thoughtful questions. In no particular order, here are my answers.
Regarding sampling rate, I did some testing with the UDB3 at rates significantly higher than 8,000 samples per cycle. I did not see any improvement in performance, so I settled at 8,000 samples per second. The design of the analog filters on the gyro outputs has an error in the UDB3 (its backwards), so there is no capacitor providing a voltage source for the sample and hold amplifiers. We have fixed on the UDB4, so maybe the UDB4 would be able to go higher. That said, I found the performance of the combination of 8,000 samples per second and "integrate and dump" to be flawless in the face of vibration. As I mentioned in my response to Andrew, MatrixPilotQuad running on a UDB3 or UDB4 has a residual in flight gyro drift rate of less than 0.5 degrees per minute.
Regarding digital gyros, I am skeptical about the way that they do filtering. Theoretical analysis of the way gyros and accelerometers are used in IMUs shows that "integrate and dump" is the right way to go for decimation filtering, and the use of a low pass filter will produce inferior results. I have not tried the digital gyros yet, but when I do I will probably want to use them to send along unfiltered samples, and then do my own filtering. Which leads into my next answer.
The decimation filter that I am using is "integrate and dump". What that means is that between executions of attitude and position computations, the outputs of the accelerometers and gyros are digitally integrated. Samples are taken at 8000 samples per second. Attitude and position are computed 40 times per second. So there are 200 samples that are integrated between computations. Each time the 40 Hz computations are executed, the digital integrators transfer their values to the computations, and then are reset. (Hence the term, "dump", in "integrate and dump".) The reason that this works so well is that the critical downstream computations are integrals. By doing the integration at 8000 samples per second, we are achieving almost the same benefits as if we were doing all of the computations at 8000 Hz. So, "integrate and dump" is a method of achieving a virtual computation rate of 8000 Hz.
Regarding the controller being active during my measurements, yes it was active. I tried turning it off, but then my quad was impossible to fly. However, the controller that I am using is not a PID. (More on that in the future.)
Regarding "impulse response", just to be clear, impulse response is a control theory term meaning the response of a system to an input that is an impulse. Once you have figured out what the impulse response of a system is, you can use it for a lot of things, including determining the Laplace transform of the system, and designing a suitable controller for the system. During my tests I never administered an impulse, but for a linear system, the total input can be thought of as a series of impulses.
Regarding the raw data, here is a link to the spreadsheet that I used to do the computations. If anyone wants to try other methods to compute the pitch impulse response of my draganflier from it have fun. Also, some of you might be interested to see how to use a spreadsheet to do the computations. Here is what is contained on the 4 sheets of the spreadsheet:
Sheet 1 : vector - matrix computations to determine impulse response
Sheet 2 : plot of h(t)
Sheet 3 : raw data
Sheet 4 : curve fitting of h(t)
Your mission, should you accept it, is to compute h(t) from pitch input and pitch output.
Pitch input is labeled pfb, for "pitch feedback". I am using a fly-by-wire system, so pfb is the only input to the differential p
-
1
-
2
-
3
of 3 Next