(or, more accurately, how we mixed our little science code into Stan and out modeling progress thus far)

This is a write-up of our experience integrating a model of linear elastic mechanical resonance in Stan.

We started off using Bayesian inference on our problem because the optimization techniques we tried weren’t working. At the recommendation of a coworker, we made a small sampler using Radford Neal’s “MCMC Using Hamiltonian Dynamics” paper (Neal 2011) and it surprised us how well it worked and how consistently reasonable the answers were. Eventually we moved to Stan to take advantage of the fancier samplers and modeling language flexibility.

We’ll save the preaching for the conclusion, but it’s our opinion that everyone out there with little science codes should be hooking them up to Stan, be they weird ordinary differential equations, complex partial differential equations, or whatever.

Formulating a science problem probabilistically usually isn’t too hard. The chief issue with getting scientific calculations into Stan are the calculations themselves. A lot of scientific codes, in one way or another, are:

  1. Too expensive to practically automatically differentiate (with the Stan autodiff)
  2. Too complicated to be of general interest to the Stan community
  3. Too technically fragile to become part of the Stan Math library

This was the case with our model, and so this notebook outlines the process of how to integrate a custom calculation (like ours) with Stan.

Shoutouts to Ben Goodrich and Bob Carpenter for walking us through this originally. The implementation here takes advantage of special interfaces in RStan and CmdStan. It’s not clear how well any of this information transfers to the other interfaces (PyStan, MatlabStan, …).

We start with a little background on our problem. It’s optional, but it should help motivate what we’re doing. We then walk through the math for a simple 1D version of our problem, how to implement this in Stan and eventually how to interface external software with Stan to solve this problem (that uses efficient, custom gradients instead of relying entirely on autodiff). Finally we go over the state of our modeling efforts and what the problems are.

Application Background

The problem we worked on was the inference bit of Resonance Ultrasound Spectroscopy (or more shortly, RUS). RUS is the process of extracting a material’s elastic properties by measuring the resonance modes of a sample of that material. how difficult an object is to stretch and deform)

Musical instruments are the standard examples of mechanical resonance at work, though the more exciting examples are bridges. The important bit is that the shape and materials of the instrument/bridge determine at what frequencies it resonates. For the most part, we know how things are shaped, and we can measure where they resonate, and so the RUS game is all about backing out the elastic constants (some standard references on this are (Visscher et al. 1991) and (Migliori et al. 1993)).

This is the Millenium Bridge in London. When it opened, feedback between peoples’ footsteps and small horizontal movements of the bridge led to visible oscillations in the structure (there are good videos of this on YouTube). Picture from Wikipedia

We don’t work with musical instruments or bridges though. The driving application here is gas turbines, as used in jet engines and land based power generators. The blades bits of the turbine that need to get hot (the blades and the rotor especially) are made of special, high temperature resistant metals called superalloys. It’s the elastic constants of superalloys that we want to know.

For superalloys, RUS gives us a few things:

  1. High precision estimates of elastic constants. Conventional mechanical testing is only accurate to 10% or so
  2. A way to evaluate the mechanical properties of materials at high operating temperatures (we’re not quite there yet, but we’re moving in this direction and are excited to see how it works)
  3. A less-destructive way to evaluate samples (you still have to destroy blades to get samples, but at least you don’t have to destroy samples to get data)

This all goes into operating turbines safely at high temperatures (which makes them more efficient).

The actual experiment works by vibrating the sample at a range frequencies and then measuring the amplitude of the response. If there is a very high peak, it is recorded as a resonance mode. This process is mostly automated. Sometimes a resonance mode does not show up, but, for the most part, the sample can be wiggled and moved around until it appears. The lowest frequency resonance modes are the hardest to measure consistently.