I came across a really interesting post yesterday about predicting when a spinning fan would come to rest. This post was very much inline with my previous Kalman filtering post (i.e. estimating and modelling a dynamic system) so I decided to give it a quick try. Instead of using the Kalman or Extended Kalman filter approach of the last post, I decided to use a least squares parameter fitting technique and try to fit a model to the data. I chose to look at two relatively simple models, one where the drag on the fan is proportional to the angular velocity and a second where the drag is proportional to the velocity squared.

In order to fit these models, I first needed to collect the tracking data. Dan provided a 30 second video of the fan spinning at what I believe is 60 frames a second. Such a long video kind of precludes manual tracking (I was not up for clicking on the screen 1500 times in a row), so I had to automate the process. As you can see in the still below, the fan blades have a few colored patches on them which help the tracking process. I chose to use the greenish colored patch closets to the light as the marker I would track.

The green marker is fixed to the blade which means the trajectory it traces through one fan revolution should maintain a constant distance from the center of the light. I used this fact to reduce the area of the image over which I tracked the object. The image below shows how I cropped the movie to help ease the tracking problem.

To actually track the green patch, I needed a way to identify green pixels and ignore non-green pixels. The normal RGB color space used in this video is not great for identifying colors. For instance, a pixel value of [0, 256, 0] is a pure green pixel, but a pixel value of [256, 256, 256] represents a pure white pixel even though the G value is identical in both pixels. A color space that allows for easier color extraction is the HSV (Hue-Saturation-Value) color space. In the HSV color space, the color (Hue) is more isolated from various lighting conditions, making color extraction more straightforward. The Hue component of the HSV representation of the above image is shown below.

The final step in identifying the green patch is determining which hue values to identify as green. This was really just a trial and error process, I tried to keep the range of Hue values as small as possible to eliminate false alarms, but large enough so that I never failed to identify the patch. In the end, I settled on using Hue values between 0.28 and 0.4 (values normalized to 1.0). The graphic below shows which pixels from the above images fall into this range.

For each green pixel, I calculated the angle about the center of the fan. Instead of taking the mean of these angles as the estimated angular displacement, I instead took the median. The rationale behind this, is that even though I tried to keep the Hue range small enough to minimize false positives, I was unable to completely eliminate them. These false positives could lead to inaccurate angular displacement estimates. If I used the mean value of the estimated angles, the false alarms would distort the estimate, by using the median angle I was able to (hopefully) exclude the one or two bad estimates. As you can see in the plot below, there do not appear to be any wildly erroneous angle estimates so the median filtering appears to have worked well.

For those that are interested, I’ve collected the displacement estimates into a text file available here.

#### Models

The equations for the models I used are shown below. The first equation represents a linear ODE where the drag force is proportional to the angular velocity. C is a coefficient representing the drag, where a larger value indicates greater drag. The second equations represents a non-linear ODE where the drag is proportional to the square of the angular velocity. Here again, C represents the drag.

The one nice thing about the system we’re modeling is that it’s free of external inputs, that is, we are observing the free response of the system. We should be able to determine a closed form solution to the free response, determine the parameters of the solution which most closely match the observed data, and then project the data into the future to determine when the fan stops. The trickiest part of the problem is determining what qualifies as a stopped fan. The equations above will never reach an angular velocity of zero by themselves. As the fan slows to a stop, it will at some point, deviate from these models and grind to a halt. I chose to classify a rotational speed of 1 revolution every 4 seconds as the speed at which the fan will most likely stop. I projected the fan’s motion into the future until the speed dropped below this threshold and marked this time as the stopping time. As a quick check of the sensitivity of the final guess to this threshold, I also looked at what time the fan would stop if the threshold had been 1 rotation every 5 seconds. I really have no idea if these thresholds are realistic but they seem reasonable.

##### Linear Model

I used an ersatz solution shown below as the general form solution to the linear drag ODE. The displacement is modeled as a decaying exponential with an initial offset.

There are three unknown parameters which must be solved for to represent this equation. I used Matlab’s lsqcurvefit algorithm to help identify these parameters. This algorithm works by searching for the parameters which will minimizes the sum of the squared data and model displacement differences. This process yielded the following values for the coefficients.

These parameters, along with the equation above, allowed me to plot the simulated displacement along with the measured displacement. I was also able to extend the data beyond the final observation to see how the model evolved.

There is a nice alignment between the data and the model within the observation period. The displacement beyond the final observed point also seems well behaved. This model can also be used to look at the angular velocity of the fan by taking the temporal derivative of the solution above. Doing this gives the graph shown below. The horizontal green line corresponds to the fan rotating at 4 seconds per turn and the blue represents a speed of 5 seconds per turn.

Using the established thresholds results in the fan coming to a stop at 38 and 43 seconds respectively.

##### Velocity Squared Model

The second model is a little more difficult than the first since it is a non-linear ODE. Luckily, the form of this non-linear ODE is well known and falls under the class of Bernoulli equations. To see this a little clearer, we can look at the original ODE in terms of the velocity instead of the displacement, this yields the following equation.

The trick to solving Bernoulli ODEs is to use a variable substitution. The following shows this substitution which yields a linear ODE that can be solved very easily.

So the ODE is now solved for in terms of the substituted variable w, we just need to use the original relationship between phi and w to get the equation in terms of phi.

One thing to note is that this phi represents velocity, not angular displacement. We need an equation in displacement to enable the data fitting. The displacement equation can be found by integrating phi with respect to time. This results in the final equation which is a solution to the original non-linear ODE above.

So here again we have three parameters that must be solved for by some means. I used the same technique as before to estimate the following parameters.

Using these parameters allowed me to again carry the fan’s motion beyond the observed data as shown below.

Qualitatively, the projected motion of the fan seems similar to the linear case, though the displacement does not seem to slow down as quickly. I sued the same velocity threshold as before to estimate the stopping time, this plot is shown below.

In this case the stopping times were approximately 61 and 79 seconds respectively. So indeed, by this model, the fan does take longer to slow down.

The graphic below compares the fit of the two methods within the data observation window. It’s difficult to make out, but the measured data is shown in black with the modeled data overlaid in dot-dashed lines. Both methods appear to fit nicely and the sum of the residuals were nearly identical, so based on this data along it’s difficult to say which method is more accurate.

#### Conclusions

Since reading the original post on Dan’s site yesterday, several commenters have posted guesses in the 70 second range. I believe the full video of the fan rotating to a stop has also been posted but I haven’t had a chance to look at it yet. It would be interesting to see how the fan behaves immediately before it stops.

As for my guess of when the fan would grind to a halt? I would tend to believe the second model where drag is proportional to velocity squared and place my bet somewhere around 70 seconds. If you look at the the displacement plot for the linear drag case, the plot seems to level off too quickly after the final observation point. There is a relatively large deceleration that takes place between 25 and 60 seconds, but this is not driven by the data, it is driven by the model. To me this seems indicative of a the model being over-constrained during the observation period leading to unrealistic behavior as we try to extrapolate the state. Take a look at the following plot of the acceleration for both models.

The acceleration of the fan changes dramatically during the first 25 seconds in the case of the second model, but in the case of the first there is not nearly as large a change. By the time the observation period has ended, the acceleration in the second model has largely leveled off, the acceleration in the first model continues to change at essentially the same rate. I don’t think this observation by itself invalidates the first model, but I am much more tempted to believe a model where the large changes happen while observing the data than a model which predicts the changes to happen after the observation period.