You want to design an experiment to measure some parameter x. Your experiment has some parameters like geometry, material properties, electronic transfer functions, etc that you can measure and control to varying degrees. Any uncertainty in your experimental parameters may lead to a systematic uncertainty in your measurement of x. We want to answer two questions:
1. Given our prior knowledge of the exerimental parameters, what is our systematic uncertainty on x?
2. What choice of experimental parameters minimizes our systematic uncertainty on x?
## Example
Currently, this repo has some code to answer (1) above.
### Experimental setup
We want to infer the mechanical loss angle of a thin film deposited on a silicon mechanical oscillator, which is clamped at its base. In this case the mechanical oscillator consists of a trident-of-tridents, which has the feature that many eigenmodes of its farthest spokes do not involve the clamped region at the base of the oscillator, protecting them from clamping losses.
We can add varying degrees of complexity to this problem, but for now assume the measurement is carried out at a fixed temperature and that the mechanical Q of the overall system is determined by a frequency-independent mechanical loss in each of the coating, substrate, and clamping material, weighted by the fraction of an eigenmode's strain-energy in each material domain.
Because each material loss contributes to eigenmode quality factor differently depending on mode shape, measuring the quality factor of many eigenmodes can tell us about the mechanical loss in the material of interest, the coating. However, our inability to perfectly predict the strain energy ratios for each eigenmode introduces a systematic error in our estimate of $`\phi_\mathrm{coating}`$.
To be concrete, our model $`A_{\vec{\mu}}`$ uses a set of parameters $`\vec{x}`$ to predict eigenfrequency and Q $`\vec{y}=\{f_i, Q_i\}`$ for the first several low-frequency eigenmodes of the mechanical system.
$`\vec{y}=A_{\vec{\mu}}(\vec{x})`$
One member of $`\vec{x}`$ is what we actually want to infer, the mechanical loss of the coating $`\phi_\mathrm{coating}`$. The other members of $`\vec{x}`$ are a variety of parameters that we may control only imperfectly (like substrate and coating thickness, or other geometric parameters), along with parameters we can only characterize imperfectly (like Young's modulus of the coating and substrate). $`A`$ may also depend on some parameters $`\vec{\mu}`$ with negligible uncertainty that we may or may not be able to control.
### Inference method
We want to use our measurements of $`\vec{y}`$ and the covariance matrix or joint probability distribution of our measurement, $`P_\mathrm{cov}(\vec{y})`$, to infer $`\vec{x}`$ with the smallest uncertainty on $`\phi_\mathrm{coating}`$ possible. We can apply Bayes' rule and the assumption that all of our model uncertainty is captured by our uncertainty in members of $`\vec{x}`$ to write
Our approach will be to use a Metropolis-Hastings algorithm to sample from $`P_\mathrm{post}`$ by drawing samples from the ratio of the known distributions $`P_\mathrm{prior}`$ and $`P_\mathrm{cov}`$.
### Simple demonstration
For
## How to run
### Guide to the files
- oneSided_Ge_coarse.m contains a COMSOL FEA model of a trident-of-tridents mechanical oscillator. The bulk of the oscillator is \<100\> crystalline silicon, it is coated with a thin film of some material of interest (in this case, I used germanium). I added a steel clamping block to the base of the oscillator, and constrained the ends of the steel block. The purpose of the model is to compute the strain energy in the clamp, silicon, and coating for each of several eigenmodes of the system.
- connect_comsol.m connects matlab to a COMSOL server session in order to run the model
- update_model.m updates the FEA model parameters
- run_experiment.m calls update_model then rebuilds the geometry and mesh, and re-solves the model
- helpers.py contains a variety of python functions needed to carry out the inference problem, including
- xinfer.py is the main script that does performs a MCMC to determine the distribution of x in the experimental parameters
### Requirements and installation
Though the python package requirements can be found in the .yml file, running COMSOL from python requires calling a MATLAB engine capable of interfacing with a COMSOL server session. Here is a description of my setup that is currently working:
1. COMSOL Multiphysics 6.0 running COMSOL Multiphysics server
2. MATLAB R2021b with python engine installed
3. Python 3.8.5
#### Installing python engine for MATLAB
The best system-dependent instructions are on [mathworks' git repo for their python engine](https://github.com/mathworks/matlab-engine-for-python/).
You can check the simplified instructions [from MathWorks](https://www.mathworks.com/help/matlab/matlab_external/install-the-matlab-engine-for-python.html), but they are somewhat misleading. You likely need to specify matlabengine version during installation, as described [on stackexchange](https://stackoverflow.com/questions/73972913/how-to-install-matlabengine-windows-10-matlab-r2020b-python-3-8-10).
### Running MCMC simulation
1. Start COMSOL server, preferably in a detachable screen or tmux session
-`comsol mphserver -login force -port 2020`
- enter a simple user and password that you don't mind appearing in plaintext. I used 'whoever' with password 'whatever'
2.`python xinfer`
- The script will print the randomly selected "true" values of the experimental parameters before running MCMC. You'll get an error if they fall outside of the prior bounds, but you may want to ensure they are satisfactory before letting your simulation run too long.
### Limitations
There are a number of complexities that should be added to the simulation. It would be beneficial to add additional geometric and material parameters to $`\vec{x}`$, such as the lengths of the tridents, Young's modulus and Poisson's ratio of the coating, temperature, and more. We should add physics-informed parametrization of mechanical loss angles, for example allowing silicon substrate mechanical loss to have a spectral or temperature dependence reflecting Akhiezer thermoelastic losses. We would also ideally perform this inference at several choices of geometric parameters to determine which ones minimize systematic uncertainty on inferred $`\phi_\mathrm{coating}`$. All of these changes expand the parameter space for MCMC, requiring more samples before the simulation converges to a reliable distribution.
Furthermore, I performed this demonstration with a coarser than optimal mesh. Preliminary mesh refinement analysis showed that using finer meshes (especially on the wide surfaces of the oscillator, and at higher eigenfrequency) does change the computed eigenfrequencies and strain energies by up to 10s%, at the expense of longer compute times.
On our local compute machine sandbox2, the COMSOL model takes \~10s to solve. Since computational power seems to be the bottleneck, the most straightforward solution is to scale up compute. emcee is an MPI-compatible python package, but we would need to install COMSOL on our hpc to run on the cluster.
