Background: The "Punkin Chunkin" contest is an annual engineering challenge in which teams use trebuchets, catapults, and air cannons in a competition to see who can hurl a pumpkin the farthest possible distance. The current record as of this writing is 4694.7 feet (1430.9 meters), nearly one mile. https://www.youtube.com/watch?v=dg4hboje3VE https://en.wikipedia.org/wiki/Pumpkin_chucking https://www.punkinchunkin.com/results/current-records.html Let's suppose we instrument a pumpkin with an accelerometer and a GPS, and we use an air cannon to launch the pumpkin straight up into the air. We can use a Kalman filter and smoother to estimate the trajectory of the pumpkin over time. For simplicity, let's suppose that the GPS and accelerometer give measurements synchronously, and that the range of the accelerometer is unlimited (although in practice there would be some upper and lower bounds, and the resulting truncation would need to be modeled). Let's consider only the vertical dimension (height). Data: The file "kf_x.txt" contains a simulated sequence of n=1990 measurements of acceleration and position (height), starting immediately after launch. The first entry in row j is the measured acceleration (meters/sec^2) and the second entry is the measured position (meters). The time between measurements is 0.01 seconds. These data were generated by simulating a true trajectory using a higher resolution physical model accounting for gravity, air resistance (using a nonlinear drag model), and random wind. The measurements were simulated as normally distributed about the true values. Model: We can use a Kalman filter and smoother to estimate the trajectory of the pumpkin. There are various ways to set up the model, but we will use what is probably the simplest possible approach. Following the notation from class and the notes, define: mu_0 = [-500 400 20] V_0 = [50^2 0 0 0 40^2 0 0 0 2^2] F = [ 1 0 0 dt 1 0 0 dt 1] H = [1 0 0 0 0 1] Q = [100*dt 0 0 0 dt 0 0 0 dt] R = [sigma_a^2 0 0 sigma_r^2] where sigma_a = 2 sigma_r = 75 dt=0.01 This choice of F encodes the physics, ignoring air resistance (velocity is the derivative of position, acceleration is the derivative of velocity), and encodes the assumption that the acceleration at the next time step should be similar to the current acceleration. This choice of H says that the measurements are unbiased. Q says that the variance of the drift/error is of order dt (which makes sense if it follows Brownian motion), and these particular values were chosen by making an initial guess for Q, doing an initial "calibration" run of the Kalman filter and smoother, and then adjusting Q to be roughly consistent with the estimated trajectory. R could be chosen using the manufacturer's specifications of the accelerometer and GPS, or by actual experimentation, or by estimating it using the calibration run; here, I just chose something reasonably close to the true values (which were sigma_a=0.98 and sigma_r=50). mu_0 is our best guess at the initial state (acceleration, velocity, and position); here, again, I just picked something reasonably close to the true values. V_0 describes our uncertainty about the initial state. Exercises: 1. Implement the Kalman filter and smoother. 2. Using the model above and the data in kf_x.txt: 2a. Run the Kalman filter, to produce mu_j for each j=1,...,n. 2b. Run the Kalman smoother, to produce mu_hat_j for each j=1,...,n. (You do not need to report your numerical results for this part.) 3. Run an IIR filter (see description below) with f=0.9 on the position data in kf_x.txt. (You do not need to report your numerical results for this part.) 4. Plot your results in the following way. Make three plots: one for acceleration, one for velocity, and one for position. On each plot, show both the output of the Kalman smoother (mu_hat_j) and the true values, which are in the file "kf_true.txt". On the position plot, also show the output of the IIR filter. On each plot, scale the time axis to be in seconds, and label each plot with the appropriate units (e.g., meters, meters/sec, etc). Discuss what you see. 5. Assess the accuracy of your results in the following way. Compute the root-mean-squared (RMS) error of each of the following, relative to the true values in kf_true.txt: 5a. acceleration, velocity, and position estimated by the Kalman filter (mu_j), 5b. acceleration, velocity, and position estimated by the Kalman smoother (mu_hat_j), 5c. position estimated by the IIR filter. So, in part 5, you will be reporting 7 numbers altogether. Include the units for each number. Discuss. 6. The model used by the Kalman filter and smoother does not account for air resistance (drag). However, it is still able to accurately estimate the trajectory (for instance, note that the pumpkin appears to reach terminal velocity on its way down). Explain. 7. Now, suppose we only had acceleration measurements (no GPS). This is called "dead reckoning". We can very easily adjust the model above to implement this by setting sigma_r to some very large value, so that the model doesn't pay any attention to the position measurements. Redo parts 4 and 5, excluding the IIR filter results, with sigma_r=10^10. Discuss. Description of IIR filter: An infinite impulse response (IIR) filter provides a simple way of taking a moving average in which more recent data points are weighted more heavily. Suppose the position measurements are r_1,...,r_n. The output of the IIR filter is the sequence s_1,...,s_n where s_1=r_1 and for j=2,...,n, s_j = f*s_{j-1} + (1-f)*r_j and f is a number between zero and one.