{ "cells": [ { "cell_type": "markdown", "id": "d81f5fb3", "metadata": {}, "source": [ "(tmcmc)=\n", "# Transitional MCMC with the classic eigenvalue problem" ] }, { "cell_type": "markdown", "id": "e0e94f8b", "metadata": {}, "source": [ "TMCMC (Transitional Markov Chain Monte Carlo) has been widely used in stochastic model updating, which aims to find the posterior samples of the parameters ($\\theta$) to be updated given empirical observations. \n", "\n", "$$\n", "\tP^j = P(\\mathcal{D}|\\boldsymbol{\\theta})^{\\beta_{j}} \\cdot P(\\boldsymbol{\\theta})\n", "$$\n", "\n", "\n", "This notebook provides a hands-on example that illustrates how *TMCMC* is used in a classic eigenvalue problem with only elementary knowledge about structural dynamics. \n", "\n", "\n", "Consider a simple two-degree-of-freedom system under undamped free vibrations. The equation of motion reads:\n", "\n", "$$\n", "\\mathbf{M} \\ddot{\\mathbf{u}}(t) + \\mathbf{K} \\mathbf{u}(t) = 0\n", "$$\n", "\n", "Assuming harmonic response yields the classic eigenvalue problem: \n", "\n", "$$\n", "\\mathbf{K} \\phi = \\lambda \\mathbf{M} \\phi\n", "$$\n", "\n", "this further leads to the characteristic equation by a nontrivial solution using the determinant:\n", "\n", "$$\n", "\\det(\\mathbf{K} - \\lambda \\mathbf{M}) = | \\mathbf{K} - \\lambda \\mathbf{M}| = 0\n", "$$\n", "\n", "where \n", "- **eigenvalues**: a system of 2 degrees of freedom will give a polynomial of 2th degree in $\\lambda$. The 2 roots will represent the frequencies at which the undamped system can oscillate in the absence of external forces, *i.e.*, natural frequencies of the system. \n", "- **eigenvectors**: The values of the vector $\\phi$ at each natural frequency describe the mode shape of the corresponding mode of vibration, which reflects the relative displacement pattern between the two masses.\n", "\n", "This setup solves the natural frequencies given known stiffness and mass matrix. Inversely, in order to identify the system given empirical measurements, consider a problem of solving the stiffness of the system given a few measurements of the $\\lambda_1$ which is contaminated by some Gaussian noise (e.g. known standard deviation $\\sigma_{\\lambda_{1}}$)." ] }, { "cell_type": "markdown", "id": "bdfee9e9", "metadata": {}, "source": [ "```{note}\n", "To further explore the identifiability of the system, Yuen (2010) proposes 3 different setups:\n", "1. Given one measured eigenvalue\n", "2. Given two measured eigenvalues\n", "3. Given a set of eigenvalue and associated eigenvector (mode shape)\n", "\n", "Readers are encouraged to have a think on the other two problems.\n", "```" ] }, { "cell_type": "markdown", "id": "4d51b48e", "metadata": {}, "source": [ "This notebook provides an exemplar solution for the problem 2 where noisy data on two eigenvalues are available to infer the stiffness parameters. Ground truth values are $\\theta_{1} = \\theta_{2} = 1.0$, and priors are assumed as $\\theta_{1} \\sim \\mathcal{U}(0.8, 2.2), \\theta_{2} \\sim \\mathcal{U}(0.4, 1.2)$." ] }, { "cell_type": "code", "execution_count": null, "id": "615fc26a-1fe6-4140-b475-58b5168ad30c", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from pyuncertainnumber.calibration import pdfs\n", "from pyuncertainnumber.calibration.tmcmc import TMCMC" ] }, { "cell_type": "code", "execution_count": null, "id": "5d533999-af6f-400a-92de-9525792956cf", "metadata": {}, "outputs": [], "source": [ "# eigen values of first mode\n", "data1 = np.array([0.3860, 0.3922, 0.4157, 0.3592, 0.3615])\n", "# eigen values of second mode\n", "data2 = np.array([2.3614, 2.5877, 2.7070, 2.3875, 2.7272])" ] }, { "cell_type": "code", "execution_count": null, "id": "e44c59f6-7c93-407c-8398-073c91d560c0", "metadata": {}, "outputs": [], "source": [ "# number of particles (to approximate the posterior)\n", "N = 1000\n", "\n", "# prior distribution of parameters\n", "k1 = pdfs.Uniform(lower=0.8, upper=2.2)\n", "k2 = pdfs.Uniform(lower=0.4, upper=1.2)\n", "\n", "# Required! a list of all parameter objects, i.e. the so-called $\\theta$ in Bayesian calibration\n", "pars_theta = [k1, k2]" ] }, { "cell_type": "markdown", "id": "ed804f8c", "metadata": {}, "source": [ "The associated likelihood is:\n", "\n", "$$\n", " p(\\mathcal{D} | \\theta) = \\frac{1}{(\\sqrt{2\\pi})^{10}\\,\\sigma_{\\lambda_1}^5\\,\\sigma_{\\lambda_2}^5}\n", "\\exp\\!\\left[\n", "-\\frac{1}{2\\sigma_{\\lambda_1}^2}\\sum_{m=1}^{5}\\left(\\tilde{\\lambda}_{1m}-\\lambda_1(\\theta)\\right)^2\n", "-\\frac{1}{2\\sigma_{\\lambda_2}^2}\\sum_{m=1}^{5}\\left(\\tilde{\\lambda}_{2m}-\\lambda_2(\\theta)\\right)^2\n", "\\right]\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "919afc11-675b-4e5c-bf0d-98a7a6150b99", "metadata": {}, "outputs": [], "source": [ "def log_likelihood_case2(particle_num, s):\n", " \"\"\"\n", " Log-likelihood for the 2DOF example - CASE 2 (two eigenvalues λ1 and λ2).\n", "\n", " Args:\n", " particle_num (int):\n", " Index of the particle (not used in this function, but required by the TMCMC framework).\n", "\n", " s (ArrayLike): numpy array of shape (2,)\n", " Parameter vector in this case:\n", " s[0] = q1 (stiffness parameter 1)\n", " s[1] = q2 (stiffness parameter 2)\n", "\n", " returns:\n", " LL (float): Log-likelihood value at this parameter vector.\n", " \"\"\"\n", " q1 = s[0]\n", " q2 = s[1]\n", "\n", " # standard deviations (5% noise on true eigenvalues λ1=0.382, λ2=2.618)\n", " sig1 = 0.0191 # 0.05 * 0.382\n", " sig2 = 0.1309 # 0.05 * 2.618\n", " center = q1 / 2.0 + q2\n", " disc = center**2 - q1 * q2\n", " if disc < 0:\n", " return -np.inf\n", "\n", " sqrt_disc = np.sqrt(disc)\n", " lambda1_s = center - sqrt_disc\n", " lambda2_s = center + sqrt_disc\n", "\n", " const_term = np.log((2 * np.pi * sig1 * sig2) ** -5)\n", " misfit1 = -0.5 * (sig1**-2) * np.sum((lambda1_s - data1) ** 2)\n", " misfit2 = -0.5 * (sig2**-2) * np.sum((lambda2_s - data2) ** 2)\n", "\n", " LL = const_term + misfit1 + misfit2\n", " return LL\n" ] }, { "cell_type": "markdown", "id": "653d9c69", "metadata": {}, "source": [ "The TMCMC algorithm is run by combining the above pieces together. " ] }, { "cell_type": "code", "execution_count": null, "id": "298bea95-0258-4f1d-911a-c5fec74ec332", "metadata": {}, "outputs": [], "source": [ "t = TMCMC(\n", " N,\n", " pars_theta,\n", " log_likelihood=log_likelihood_case2,\n", " status_file_name=\"status_file_2DOF_problem2.txt\",\n", ")" ] }, { "cell_type": "markdown", "id": "7d63b709", "metadata": {}, "source": [ "```{warning}\n", "Do not run the code below in a Jupyter environment !!! Instead, run the whole code as a Python file. This is due to the incompatbility of the \"multiprocess\" package in the Jupyter environment. This noteboos is meant as a demonstration of how the code can be used. Unfortunately, due to the incompatbility, the code cannot directly run with the notebook. Readers are encouraged to copy all the code into aother python file to test the code. \n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "78d98520", "metadata": {}, "outputs": [], "source": [ "mytrace = t.run() # Do not run in a Jupyter environment due to multiprocessing issues" ] }, { "cell_type": "markdown", "id": "acce6f42", "metadata": {}, "source": [ "
\n", " \n", "
Posterior of the problem 2
\n", "
\n", "\n", "A trace object will be returned. Meanwhile, a log file named \"\"status_file_2DOF_problem2.txt\" containing information about the sampling process will be saved." ] } ], "metadata": { "kernelspec": { "display_name": "pbox_nn", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.10" } }, "nbformat": 4, "nbformat_minor": 5 }