Case Study

Cantilever beam

Thyge Vinther Ludvigsen, s203591

29 Apr, 2026

Time Plan

New Stochastic model

Sensor error model

Multivariate normal distribution:

\[ \epsilon_{\text{OMA}} \sim \mathcal{N}(0, \boldsymbol\Sigma_{OMA}(e)) \]

Diagonal covariance matrix: 1

\[ \Sigma_{\text{OMA}}(e)= \begin{bmatrix} \Sigma_{\omega}(e) & 0 \\ 0 & \Sigma_{\phi}(e) \end{bmatrix} \]

Natural frequencies

Frequency covariance Matrix: \[ \Sigma_{\omega}(e) = \operatorname{diag}\!\left( \operatorname{Var}(\hat{\omega}_1 \mid e), \dots, \operatorname{Var}(\hat{\omega}_m \mid e) \right) \]

Error approximation for frequencies:

\[ \boxed{ \operatorname{Var}(\hat{\omega}_i \mid e) \;\approx\; \frac{\sigma_\varepsilon^2} {\;\|\mathbf S(e) \odot \phi_i\|_2^2} } \]

Where:

  • \(\sigma_\varepsilon = CV \cdot \mu_{\omega_i}\): sensor noise, given as sensor CV (coefficient of variation) multiplied by the mean value of the natural frequencies.
  • \(\|\mathbf S(e) \odot \phi_i\|_2^2\): is the modal observability,
    • \(\| x \|_2\) is the L2 norm of the vector \(x\), so the euclidean distance.
    • \(\mathbf S(e)\) is the sensor configuration vector.
    • \(\phi_i\) is the mode shape vector for the \(i\)-th mode, which describes the deformation pattern of the structure at that mode.
    • \(\odot\) denotes the element-wise (Hadamard) product, which multiplies corresponding elements of the two vectors.

\[ \|\mathbf S(e) \odot \phi_i\|_2 = \sqrt{\sum_{j=1}^{n} (S_j(e) \cdot \phi_{i,j})^2} \]

\[ \sum_{j=1}^{n} (S_j(e) \cdot \phi_{i,j})^2 = (S_1(e) \cdot \phi_{i,1})^2 + (S_2(e) \cdot \phi_{i,2})^2 + ... + (S_n(e) \cdot \phi_{i,n})^2 \]

S = np.array([0, 1, 1])
phi = np.array([-1, 0.2, 0.5])
print(f"np.multiply(S, phi) = {np.multiply(S, phi)}")
print(f"np.linalg.norm(np.multiply(S, phi)) = {np.linalg.norm(np.multiply(S, phi))}") # element-wise multiplication and then take the L2 norm
print(f"Manual calculation: {np.sqrt(0.5**2 + 0.2**2)}") # manual calculation of the L2 norm
np.multiply(S, phi) = [-0.   0.2  0.5]
np.linalg.norm(np.multiply(S, phi)) = 0.5385164807134505
Manual calculation: 0.5385164807134505

Mean values of |phi|: 0.4937138672223707
Number of samples: 1000

Mean values of |phi|: 0.49910643511437275
Number of samples: 1000

Mean values of |phi|: 0.5022195726544022
Number of samples: 1000

Mean values of |phi|: 0.503364421468163
Number of samples: 1000

Mode shapes

Mode shape covariance matrix: 1

\[ \Sigma_{\phi}(e) = \operatorname{diag}\!\left( \operatorname{Var}(\hat{\phi}_1 \mid e), \dots, \operatorname{Var}(\hat{\phi}_m \mid e) \right) \]

Error approximation for mode shapes: \[ \boxed{ \operatorname{Var}(\hat{\phi}_i \mid e) \;\approx\; \frac{ \sigma_\varepsilon^2 } {\;\|\mathbf S(e) \odot \phi_i\|_2^{}} } \]

Where:

  • \(\sigma_\varepsilon = \sigma_\phi\): where \(\sigma_\phi\) is the sensor noise for mode shapes I set it to be the same for all modes.
  • \(\|\mathbf S(e) \odot \phi_i\|_2\): is the modal observability.
    • \(\| x \|_2\) is the L2 norm of the vector \(x\), so the euclidean distance.
    • \(\mathbf S(e)\) is the sensor configuration vector.
    • \(\phi_i\) is the mode shape vector for the \(i\)-th mode, which describes the deformation pattern of the structure at that mode.
    • \(\odot\) denotes the element-wise (Hadamard) product, which multiplies corresponding elements of the two vectors.

Features

Feature vector:

\[ \boxed{ \mathbf y = h(\bar{\boldsymbol\phi}, \bar{\boldsymbol\omega}) } \]

Where \(h(\cdot)\) is an function that maps the modal parameters to the feature space.

  • Natural frequencies, \(\omega\)
  • Total Modal Assurance Criterion (TMAC)
    • \(MAC(\phi_a,\phi_b)=\frac{|\phi_a^T\phi_b|^2}{(\phi_a^T\phi_a)(\phi_b^T\phi_b)}\)
    • \(TMAC(H_{\text{no damage}}, \bar{\boldsymbol\phi}) = \sum_{i=1}^{n} \frac{|\phi_{0,i}^T \bar{\phi}_i|^2}{(\phi_{0,i}^T \bar{\phi}_i)^2}\)
  • Modal Flexibility: \(F = \sum_{i=1}^{n} \frac{\phi_i \phi_i^T}{\omega_i^2}\)

Problems

  • Features with nonlinear combinations of stochastic variables leads to implicit likelihoods.

Code example - Stochastic model

Goal of this code example

  1. Setup an new stochastic model
  2. generate modal parameters given an sensor configuration.
  3. Compute the features for an given sensor configuration.

Give:

  • Modal parameters for all states, \(\omega\) and \(\phi\)
  • sensor coefficients of variation (CV)
  • Sensor configuration vector, \(S\)
  • Prior probabilities for the states, \(p(\theta)\)

Test:

  • Calculate the variance of the modal parameters given \(S\) and CV
  • Generate samples of modal parameters given \(S\) and CV
  • Compute the features for a given sensor configuration.
import numpy as np
from numpy.random import multivariate_normal

class BR:
  def __init__(self, omega_all, phi_all, CV_omega, sigma_phi, prior, N = 10**4, print_code = False):
    # all parameters for later use
    self.omega_all = omega_all
    self.phi_all = phi_all
    # No damage parameters
    self.omega_H0 = omega_all[0]
    self.phi_H0 = phi_all[0]
    # Damage parameters
    self.omega_H1 = omega_all[1:]
    self.phi_H1 = phi_all[1:] 
    
    # Coefficient of variation for omega and phi
    self.CV_omega = CV_omega
    self.sigma_phi = sigma_phi

    # number of Monte Carlo samples
    self.N = N
    
    # prior probabilities for theta
    self.prior = prior

    # print code
    self.print_code = print_code


  # def feature_vector(self, features_vector):

  #   return 

  # Statical parameters for the multivariate normal distribution 
  def covariance_matrix(self, S):
    # Sensor vector
    S = np.array(S)

    # Standard deviation
    sigma_omega = self.CV_omega * self.omega_H0 
    sigma_phi = self.sigma_phi * np.ones(self.phi_H0.shape[0]) 

    # modal observability
    mo = np.array([np.linalg.norm(np.multiply(S, phi_i)) for phi_i in self.phi_H0])

    # Variance Omega
    var_omega = sigma_omega**2 / (mo**2)

    # Variance Phi
    var_phi = sigma_phi**2 / (mo**2)

    # Diagonal covariance matrix
    # mu = np.concatenate([self.omega_H0, self.mask_and_normalize(self.phi_H0, S).flatten()])  # mean vector
    mu = self.mu_vector(self.omega_H0, self.phi_H0, S) # mean vector
    diagonal = np.concatenate([var_omega, np.repeat(var_phi, np.sum(S, dtype=int))])  # variances for omega and phi
    cov = np.diag(diagonal)  # covariance matrix
    
    # save for later use
    self.mu = mu
    self.cov = cov 
  
  def mu_vector(self, omega, phi, S):
    # Sensor vector
    S = np.array(S)
    # mean vector 
    mu = np.concatenate([omega, self.mask_and_normalize(phi, S).flatten()])  # mean vector
  
    return mu

  def mask_and_normalize(self,X, S):
    mask = S == 1

    # Step 1: apply mask on last axis
    filtered = X[..., mask]

    # Step 2: compute max along last axis
    max_vals = np.max(np.abs(filtered), axis=-1, keepdims=True)

    # Step 3: avoid division by zero
    max_vals[max_vals == 0] = 1

    # Step 4: normalize
    return filtered / max_vals


  def MCS_modal(self, S, prior_Importance = None):
    # S is the sensor vector
    S = np.array(S)

    # get covariance matrix and mean vector
    self.covariance_matrix(S)

    if prior_Importance is not None:
      # Imoportance sampling for structural states: theta ~ Q(theta) where Q is the importance distribution
      states = np.random.choice(range(len(self.prior)), size=self.N, p=prior_Importance)
      H_vec = (states > 0).astype(int)  # 0 if state == 0, else 1
      self.states = states
      self.H_vec = H_vec

    elif prior_Importance is None:
      # Generate all states and H labels at once: theta ~ P(theta)
      states = np.random.choice(range(len(self.prior)), size=self.N, p=self.prior)
      H_vec = (states > 0).astype(int)  # 0 if state == 0, else 1
      self.states = states
      self.H_vec = H_vec

    # generate errors for omega and phi
    # errors = np.random.multivariate_normal(mean=np.zeros(self.cov.shape[0]), cov=self.cov, size=self.N)

    # model properties 
    # lambda_star = np.random.multivariate_normal(mean=self.mu, cov=self.cov, size=self.N)
    lambda_star = np.random.multivariate_normal(mean=np.zeros(self.cov.shape[0]), cov=self.cov, size=self.N)
    lambda_star[:, :len(self.omega_H0)] += self.omega_all[states]
    lambda_star[:, len(self.omega_H0):] += self.mask_and_normalize(self.phi_all[states], S).reshape(self.N, -1)
    self.lambda_star = lambda_star

    # extract omega and phi from lambda_star
    omega_bar = lambda_star[:, :self.omega_H0.shape[0]]
    # phi_bar = lambda_star[:, self.omega_H0.shape[0]:]
    phi_bar = lambda_star[:, self.omega_H0.shape[0]:].reshape(self.N, self.phi_H0.shape[0], np.sum(S, dtype=int))
    self.phi_no_norm = phi_bar.copy() # save the non-normalized phi for later use
    # normize phi_bar by  by the maximum absolute value across all samples for each mode
    # normalize each mode-shape vector per sample (along the last axis)
    max_abs = np.max(np.abs(phi_bar), axis=2, keepdims=True)  # shape: (N, N_phi, 1)
    phi_bar = phi_bar / np.where(max_abs == 0, 1, max_abs)
    # phi_bar = self.mask_and_normalize(phi_bar, S)

    self.omega_bar = omega_bar
    self.phi_bar = phi_bar
    

  def MAC(self, phi1, phi2):
    # Compute the Modal Assurance Criterion (MAC) between two mode shapes
    numerator = np.abs(np.dot(phi1, phi2))**2
    denominator = np.dot(phi1, phi1) * np.dot(phi2, phi2)
    return numerator / denominator if denominator != 0 else 0

  def TMAC(self, phi_theta_1, phi_theta_2):
    N_samples = phi_theta_1.shape[0]
    TMAC_matrix = np.zeros((N_samples, N_samples))
    for i in range(N_samples):
        for j in range(N_samples):
            TMAC_matrix[i, j] = self.MAC(phi_theta_1[i], phi_theta_2[j])

    # sum of diagonal of the MAC 
    self.print2("TMAC matrix:\n", TMAC_matrix)  
    TMAC = np.sum(np.diag(TMAC_matrix)) / (N_samples)
    return TMAC
  
  # Modal Flexibility (MF)
  def MF(self, phi_vec, omega_vec):
    # Compute the Modal Flexibility (MF) for a given mode shape and natural frequency
    return sum(np.diag(np.abs(np.dot(phi_vec, phi_vec.T) / (omega_vec**2))))

  def print2(self, *args):
    if self.print_code:
        print(*args)

Code - Setup

np.random.seed(32)

# frequencyes and mode shapes
omega_all = np.array([[1, 6.72223044],    # H0, theta_0
                      [1.2, 8.72223044],  # H1, theta_1
                      [1.6, 10.8]])          # H1, theta_2
phi_all = np.array([[[0.4, 0.5, 1], [-0.3, -1, .5]],   # H0, theta_0
                    [[0.2, 0.2, 1], [-0.5, -1, .5]],   # H1, theta_1
                    [[-0.1, 0.1, 1], [-0.9, -1, .5]]]) # H1, theta_2

# Coefficient of variation (CV)
CV_omega = 0.1
sigma_phi = 0.2

# Prior probabilities for the states (theta)
prior = np.array([0.5, 0.4, 0.1])

# Number of Samples 
N = 10**4

# Instantiate the BR class
BRinf = BR(omega_all, phi_all, CV_omega, sigma_phi, prior, N)

# Sensor configuration
S = np.array([1, 1, 1])

# Generate modal parameters given the sensor configuration
BRinf.MCS_modal(S)

"""Simple results"""
print("Sampled states (theta):", BRinf.states)
print(f"#θ_0 = {sum(BRinf.states == 0)}, fraction = {sum(BRinf.states == 0) / BRinf.N:.2f}")
print(f"#θ_1 = {sum(BRinf.states == 1)}, fraction = {sum(BRinf.states == 1) / BRinf.N:.2f}")
print(f"#θ_2 = {sum(BRinf.states == 2)}, fraction = {sum(BRinf.states == 2) / BRinf.N:.2f}")

print("\nSampled modal parameters")
print("-" * 40)
print("Frequencies:")
print(np.round(BRinf.omega_bar[0], 4))

print("\nMode shapes:")
print(np.round(BRinf.phi_bar[0], 4))
Sampled states (theta): [1 0 1 ... 1 0 1]
#θ_0 = 4915, fraction = 0.49
#θ_1 = 4071, fraction = 0.41
#θ_2 = 1014, fraction = 0.10

Sampled modal parameters
----------------------------------------
Frequencies:
[1.2598 8.996 ]

Mode shapes:
[[ 0.3688  0.1529  1.    ]
 [-0.4113 -1.      0.3475]]

Plot - Omega

Plot - Mode shapes 1

Plot - MAC + TMAC

Plot - Modal Flexibility (MF)

case study Overview

Goal of this Cass study

  • Setup an new stochastic model
  • Evaluate the Bayes risk whit the new stochastic model
  • Try out new decision rules
  • Implement importance sampling in MCS

Chosen Features

Modal properties:

  • x first natural frequencies
  • x first mode shapes

Features:

  • natural frequencies
  • Total Modal Assurance Criterion (TMAC)
    • \(MAC(\phi_a,\phi_b)=\frac{|\phi_a^T\phi_b|^2}{(\phi_a^T\phi_a)(\phi_b^T\phi_b)}\)
    • with no damage as reference \(H_0\)
  • Modal Flexibility
    • \(F = \sum_{i=1}^{n} \frac{\phi_i \phi_i^T}{\omega_i^2}\)

Evaluation of the Bayes risk

Numarical integration:

  • Monte Carlo simulation (MCS)
  • MCS with importance sampling

Decision rules:

  • Explicit likelihood for Modal properties (no use of features)
    • Maximum a posteriori (MAP) decision rule \(\delta(y) = \arg\max_{a} p(H \mid y)\)
  • Supervised Machine learning (ML) classification based on the features (no use of likelihood)
    • Logistic regression
    • Tree-based methods

Prior distibution

H_1 = 0.01
H_0 = 1-H_1

# uninformed prior
N_prior = 2**4
prior_H = np.array([H_0, H_1])
prior_theta = np.concatenate(([H_0], np.repeat(H_1 / (N_prior - 1), N_prior - 1)))

print(f"Prior probabilities for H: {prior_H}")
print(f"\nPrior probabilities for θ: {np.round(prior_theta, 4)}")
Prior probabilities for H: [0.99 0.01]

Prior probabilities for θ: [9.9e-01 7.0e-04 7.0e-04 7.0e-04 7.0e-04 7.0e-04 7.0e-04 7.0e-04 7.0e-04
 7.0e-04 7.0e-04 7.0e-04 7.0e-04 7.0e-04 7.0e-04 7.0e-04]

Cost function

\[ C_\text{total} = C_\text{prediction}(H, d) + C_\text{sensors}(e) \]

\[ C_\text{sensors}(e) = 5 + 1 \cdot \text{number of sensors} \]

\[ C_\text{prediction} = \begin{bmatrix} 0 & 50 \\ 10000 & 10 \end{bmatrix} = \begin{bmatrix} \text{True Negative} & \text{False Positive} \\ \text{False Negative} & \text{True Positive} \end{bmatrix} \]

where:

  • Positive: is the prediction of damage.
  • Negative: is the prediction of no damage.
  • True Positive (TP): Correctly predicting damage when it is present.
  • False Positive (FP): Incorrectly predicting damage when it is not present.
  • True Negative (TN): Correctly predicting no damage when it is not present.
  • False Negative (FN): Incorrectly predicting no damage when there is actually damage.
# C [H, d] = the prediction cost
prediction_cost = np.array([[0,   50],  # TN, FP
                            [10000, 100]   # FN, TP
                            ])  

H = 1; d = 1 
print(f"Prediction cost: {prediction_cost[H, d]}, which is the cost of predicting {d} when the true state is {H}.")

def sensor_cost(n_sensors):
    return 5 + 1 * n_sensors
Prediction cost: 100, which is the cost of predicting 1 when the true state is 1.

FEM model

Geometry

Mode Shapes

Damaged scenarios

Damage scenarios: 
[(0, 0, 0, 0), (0, 0, 0, 1), (0, 0, 1, 0), (0, 0, 1, 1), (0, 1, 0, 0), (0, 1, 0, 1), (0, 1, 1, 0), (0, 1, 1, 1), (1, 0, 0, 0), (1, 0, 0, 1), (1, 0, 1, 0), (1, 0, 1, 1), (1, 1, 0, 0), (1, 1, 0, 1), (1, 1, 1, 0), (1, 1, 1, 1)]
size of theta: 16
Plot of case: (1, 1, 1, 1):

Sensor configurations

# nodes  
e_sensor = np.arange(len(X))

# list of fixed nodes
fixed_nodes = [bounds[i][0] for i in range(len(bounds))]
fixed_nodes
# remove fixed nodes
e_sensor = e_sensor[~np.isin(e_sensor, fixed_nodes)] # ~ is the logical NOT operator, which inverts the boolean array returned by np.isin
print(f"Sensor nodes: {e_sensor}")

# make an list whit all possible combinations of sensor placements
import itertools # generate the solution space of all possible combinations.
# e_configuration = list(itertools.combinations(e_sensor.tolist(), 3))
e_configuration =  np.array(list(itertools.product([0, 1], repeat=len(e_sensor)))) # all combinations of 0 and 1 for each sensor node
e_configuration = e_configuration 
print(f"Number of sensor configurations: {len(e_configuration)}")
for i in range(len(e_configuration)):
    print(f"Configuration {i} given in nodes with sensors: {e_configuration[i]}")
Sensor nodes: [1 2 3 4]
Number of sensor configurations: 16
Configuration 0 given in nodes with sensors: [0 0 0 0]
Configuration 1 given in nodes with sensors: [0 0 0 1]
Configuration 2 given in nodes with sensors: [0 0 1 0]
Configuration 3 given in nodes with sensors: [0 0 1 1]
Configuration 4 given in nodes with sensors: [0 1 0 0]
Configuration 5 given in nodes with sensors: [0 1 0 1]
Configuration 6 given in nodes with sensors: [0 1 1 0]
Configuration 7 given in nodes with sensors: [0 1 1 1]
Configuration 8 given in nodes with sensors: [1 0 0 0]
Configuration 9 given in nodes with sensors: [1 0 0 1]
Configuration 10 given in nodes with sensors: [1 0 1 0]
Configuration 11 given in nodes with sensors: [1 0 1 1]
Configuration 12 given in nodes with sensors: [1 1 0 0]
Configuration 13 given in nodes with sensors: [1 1 0 1]
Configuration 14 given in nodes with sensors: [1 1 1 0]
Configuration 15 given in nodes with sensors: [1 1 1 1]

Covariance matrix

Plot’s

Covariance matrix:
 [1.000e-02 5.500e-01 3.200e+00 2.842e+01 1.100e-01 1.100e-01 1.900e-01
 1.900e-01 1.400e-01 1.400e-01 3.900e-01 3.900e-01]
Sensor configuration: [0 0 1 1]

Covariance matrix:
 [0.01 0.22 1.61 4.98 0.1  0.1  0.1  0.1  0.08 0.08 0.08 0.08 0.07 0.07
 0.07 0.07 0.07 0.07 0.07 0.07]
Sensor configuration: [1 1 1 1]

Analysis 1

  • Evaluation method: Monte Carlo simulation (MCS)
  • likelihood: Explicit likelihood for Modal properties (no use of features)
  • Decision rule: Maximum a posteriori (MAP) decision rule \(d(y) = \arg\max_{H} p(H \mid y)\)

Workflow - Normal

%%{init: {'flowchart': {'useMaxWidth': true}}}%%

flowchart TB
    subgraph SensorDesign["Sensor configurations: e    ."]
        direction TB
        style SensorDesign fill:#e8f5e9,stroke:#2e7d32,stroke-width:3px,color:#1b5e20
        
            
        subgraph Measurement["Monte Carlo Simulations: N    ......"]
            direction TB
            style Measurement fill:#fff3e0,stroke:#e65100,stroke-width:3px,color:#bf360c

            D2 --> B

            A2{"Structural strategy (MCS)<br/> θ* ~ P(θ)"}
            B2{"Generate error (MCS)<br/> ε(e) ~ f(ε | e)"}
            C2["Features<br/>y = g(θ*, e) + ε(e)"] 

            

            A2 --> B2
            B2 --> C2
            C2 --> D2

            style C2 D2 fill:#ffebee,stroke:#c62828,stroke-width:2px,color:#b71c1c
            style D2 fill:#ffebee,stroke:#c62828,stroke-width:2px,color:#b71c1c
            style B fill:#ffebee,stroke:#c62828,stroke-width:2px,color:#b71c1c

            subgraph StructuralState["Structural states: θ    ."]
                direction TB
                style StructuralState fill:#e3f2fd,stroke:#1565c0,stroke-width:3px,color:#0d47a1

                D2{"Likelihood of observation (MCS)<br/>f(y | θ)"}
                B["Likelihood<br/>f(y | H) = ∑<sub>θ</sub> f(y | θ) P(θ | H)"]
                end

        subgraph DamageNoDamage["Structural categories: H    ."]
            direction TB
            style DamageNoDamage fill:#fff9c4,stroke:#f57f17,stroke-width:3px,color:#e65100
        
            C["Posterior<br/>P(H | y)"]
            D["Decision Rule<br/>d = max P(H | y)"]
            end
        E["Loss function<br/> L(H, d) "]

        B --> C
        C --> D
        D --> E
        F["Bayes Risk<br/> Ψ(e) =1/N  ∑<sub>N</sub>  L(H, d)"]


        end

        E --> F


        H["Optimization<br/>e* = argmin Ψ(e)"]
        F --> H
    end

Workflow - Explicit likelihood

%%{init: {'flowchart': {'useMaxWidth': true}}}%%
flowchart TB
    subgraph SensorDesign["Sensor configurations: e    ."]
        direction TB
        style SensorDesign fill:#e8f5e9,stroke:#2e7d32,stroke-width:3px,color:#1b5e20
        
            
        subgraph Measurement["Monte Carlo Simulations: N    ......"]
            direction TB
            style Measurement fill:#fff3e0,stroke:#e65100,stroke-width:3px,color:#bf360c

            D2 --> B

            A2{"Structural strategy (MCS)<br/> θ* ~ P(θ)"}
            B2{"Generate error (MCS)<br/> ε(e) ~ f(ε | e)"}

            

            A2 --> B2
            B2 --> D2

            style D2 fill:#ffebee,stroke:#c62828,stroke-width:2px,color:#b71c1c
            style B fill:#ffebee,stroke:#c62828,stroke-width:2px,color:#b71c1c

            subgraph StructuralState["Structural states: θ    ."]
                direction TB
                style StructuralState fill:#e3f2fd,stroke:#1565c0,stroke-width:3px,color:#0d47a1

                D2{"Likelihood of observation (MCS)<br/>f(ε(e) | θ)"}
                B["Likelihood<br/>f(ε(e) | H) = ∑<sub>θ</sub> f(ε(e) | θ) P(θ | H)"]
                end

        subgraph DamageNoDamage["Structural categories: H    ."]
            direction TB
            style DamageNoDamage fill:#fff9c4,stroke:#f57f17,stroke-width:3px,color:#e65100
        
            C["Posterior<br/>P(H | ε(e))"]
            D["Decision Rule<br/>d = max P(H | ε(e))"]
            end
        E["Loss function<br/> L(H, d) "]

        B --> C
        C --> D
        D --> E
        F["Bayes Risk<br/> Ψ(e) =1/N  ∑<sub>N</sub>  L(H, d)"]


        end

        E --> F


        H["Optimization<br/>e* = argmin Ψ(e)"]
        F --> H
    end

Code - Setup

from scipy.stats import multivariate_normal as mvn
from sklearn.metrics import confusion_matrix

# frequencyes and mode shapes
omega_all = load["natural_frequencies"]
phi_all = mode_shapes

# Coefficient of variation (CV)
CV_omega = 0.1
sigma_phi = 0.1

# Prior probabilities for the states (theta)
prior_theta = prior_theta

# Number of Samples 
N = 10**4

# must have an minimum of 2 sensors to be able to detect damage
e_config_2 = e_configuration[e_configuration.sum(axis=1) >= 2]


theta = OMA.theta
bayes_risk_sensors = np.zeros(len(e_config_2))
cm = np.zeros((len(e_config_2), 2, 2), dtype=int) # confusion matrix for each sensor configuration

# Instantiate the BR class
A1 = BR(omega_all, phi_all, CV_omega, sigma_phi, prior_theta, N)

# MCS for sensor configuration S
# S = e_configuration[-1]
# print(f"Running MCS for sensor configuration: {S}")

# print(f"e_configuration[1:]:\n {e_configuration[1:]}")
for ei, S in enumerate(e_config_2):
    # Generate samples
    A1.MCS_modal(S)

    # print(f"Running MCS for sensor configuration: {S}")

    # Vectorized Likelihood calculation:  pdf = f(y|theta_i)
    pdf_epsilon_all = np.zeros((N, len(theta)))
    for i in range(len(theta)):
        mu = A1.mu_vector(omega_all[i], phi_all[i], S)
        pdf_epsilon_all[:, i] = mvn.pdf(A1.lambda_star - mu, mean=np.zeros_like(mu), cov=A1.cov)

    # Likelihood: f(y|H) = sum_i f(y|theta_i) * P(theta_i|H)
    Likelihood_H0 = pdf_epsilon_all[:, 0] * prior_theta[0] / prior_H[0] # P(theta_i|H) = prior[0] / P_H0 for the undamaged state
    # Sum across all damage states (i=1 to end)
    Likelihood_H1 = np.sum(pdf_epsilon_all[:, 1:] * (prior_theta[1:] / prior_H[1]), axis=1)

    # Posterior probabilities: P(H|y) ∝ f(y|H) * P(H)
    Post0 = Likelihood_H0 * prior_H[0]
    Post1 = Likelihood_H1 * prior_H[1]

    # Decisions rule: maximimum a posteriori (MAP) 
    d_vec = (Post1 > Post0).astype(int)

    # Loss and Bayes Risk
    total_costs = prediction_cost[A1.H_vec, d_vec] + sensor_cost(len(e_configuration[sum(S)]))

    # Confusion matrix
    cm[ei] = confusion_matrix(A1.H_vec, d_vec)

    # Estimation of the bayes risk using Monte Carlo.
    Bayes_Risk = np.mean(total_costs)

    print(f"Bayes Risk: {Bayes_Risk}")
    bayes_risk_sensors[ei] = Bayes_Risk

Bayes risk for different sensor configurations

Minimum Bayes Risk: 14.095
Best sensor configuration: [1 1 1 0] (index: 9)

Confusion matrices - Score’s

Confusion matrices: \[ \text{Confusion Matrix} = \begin{bmatrix} TN & FP \\ FN & TP \end{bmatrix} \]

Recall: \[ recall = \frac{TP}{TP + FN} \]

Precision: \[ precision = \frac{TP}{TP + FP} \]

F1-score: \[ F1 = 2 \cdot \frac{\text{precision} \cdot \text{recall}}{\text{precision} + \text{recall}} \]

Confusion matrices - Plot 1

Confusion matrices - Plot 2

Note

There is an need for Importance sampling

Down sides

Bayesian inference whit the likelihood function:

  1. I assume perfect knowledge of the statistical parameters for the error model for an given sensor configuration?

Analysis 1.1

  • Evaluation method: Monte Carlo simulation with importance sampling
  • likelihood: Explicit likelihood for Modal properties (no use of features)
  • Decision rule: Maximum a posteriori (MAP) decision rule \(d(y) = \arg\max_{H} p(H \mid y)\)

Importance sampling

Variance reduction of expectation: 1

\[ E[g(X)] = \int g(x) f(x) dx = \int g(x) w(x) h(x) dx = E_h[g(X) w(X)] \simeq \frac{1}{N} \sum_{i=1}^{N} g(x_i) w(x_i) \]

Where:

  • \(f(x)\) is the probability density function (PDF) of the original distribution from which we want to sample.
  • \(g(x)\) is the function we want to compute the expectation of.
  • \(h(x)\) is the PDF of the proposal distribution.
  • \(w(x) = \frac{f(x)}{h(x)}\) is the importance weight, which corrects for the fact that we are sampling from \(h(x)\) instead of \(f(x)\).

Importance sampling code

from scipy.stats import multivariate_normal as mvn
from sklearn.metrics import confusion_matrix

# frequencyes and mode shapes
omega_all = load["natural_frequencies"]
phi_all = mode_shapes

# Coefficient of variation (CV)
CV_omega = 0.1
sigma_phi = 0.1

# Prior probabilities for the states (theta)
prior_theta = prior_theta

# Number of Samples 
N = 10**4

# must have an minimum of 2 sensors to be able to detect damage
e_config_2 = e_configuration[e_configuration.sum(axis=1) >= 2]


theta = OMA.theta
bayes_risk_sensors = np.zeros(len(e_config_2))
cm = np.zeros((len(e_config_2), 2, 2), dtype=int) # confusion matrix for each sensor configuration

# Instantiate the BR class
A1 = BR(omega_all, phi_all, CV_omega, sigma_phi, prior_theta, N)

# improtance sampling prior for the states (theta) H_0 = 0.5 and H_1 = 0.5
prior_H_Importance = np.array([0.5, 0.5])
prior_Importance = np.concatenate(([prior_H_Importance[0]], np.repeat(prior_H_Importance[1] / (len(prior_theta) - 1), len(prior_theta) - 1)))

print(f"Importance sampling prior probabilities for H: {prior_H_Importance}")
print(f"\nImportance sampling prior probabilities for θ: {np.round(prior_Importance, 4)}")

# print(f"e_configuration[1:]:\n {e_configuration[1:]}")
for ei, S in enumerate(e_config_2):
    # Generate samples
    A1.MCS_modal(S, prior_Importance=prior_Importance)

    # print(f"Running MCS for sensor configuration: {S}")

    # Vectorized Likelihood calculation:  pdf = f(y|theta_i)
    pdf_epsilon_all = np.zeros((N, len(theta)))
    for i in range(len(theta)):
        mu = A1.mu_vector(omega_all[i], phi_all[i], S)
        pdf_epsilon_all[:, i] = mvn.pdf(A1.lambda_star - mu, mean=np.zeros_like(mu), cov=A1.cov)

    # Likelihood: f(y|H) = sum_i f(y|theta_i) * P(theta_i|H)
    Likelihood_H0 = pdf_epsilon_all[:, 0] * prior_theta[0] / prior_H[0] # P(theta_i|H) = prior[0] / P_H0 for the undamaged state
    # Sum across all damage states (i=1 to end)
    Likelihood_H1 = np.sum(pdf_epsilon_all[:, 1:] * (prior_theta[1:] / prior_H[1]), axis=1)

    # Posterior probabilities: P(H|y) ∝ f(y|H) * P(H)
    Post0 = Likelihood_H0 * prior_H[0]
    Post1 = Likelihood_H1 * prior_H[1]

    # Decisions rule: maximimum a posteriori (MAP) 
    d_vec = (Post1 > Post0).astype(int)

    # Loss and Bayes Risk
    total_costs = prediction_cost[A1.H_vec, d_vec] + sensor_cost(len(e_configuration[sum(S)]))

    # Confusion matrix
    cm[ei] = confusion_matrix(A1.H_vec, d_vec)

    # Importance weights for the samples
    w = prior_H/prior_H_Importance
    importance_weights = w[A1.H_vec]  # Get the importance weights for each sample based on its true state H_vec


    # Estimation of the bayes risk using Monte Carlo.
    Bayes_Risk = np.mean(importance_weights * total_costs)

    print(f"Bayes Risk: {Bayes_Risk}")
    bayes_risk_sensors[ei] = Bayes_Risk
Importance sampling prior probabilities for H: [0.5 0.5]

Importance sampling prior probabilities for θ: [0.5    0.0333 0.0333 0.0333 0.0333 0.0333 0.0333 0.0333 0.0333 0.0333
 0.0333 0.0333 0.0333 0.0333 0.0333 0.0333]
Bayes Risk: 52.273048
Bayes Risk: 47.570064
Bayes Risk: 43.8283
Bayes Risk: 32.355183999999994
Bayes Risk: 50.633292000000004
Bayes Risk: 35.66810400000001
Bayes Risk: 29.192908000000003
Bayes Risk: 36.016716
Bayes Risk: 27.297348
Bayes Risk: 20.894144
Bayes Risk: 15.547128

Bayes risk for different sensor configurations

Note

Compare this with the results from non-importance sampling

Minimum Bayes Risk: 15.547128
Best sensor configuration: [1 1 1 1] (index: 10)

Confusion matrices - Plot 1

Note

Compare this with the results from non-importance sampling

Confusion matrices - Plot 2

Note

Now is it good

Analysis 2

  • Evaluation method: Monte Carlo simulation (MCS)
  • likelihood: No likelihood, use supervised machine learning (ML) classification based on the features.
  • Decision rule: use the trained ML model to directly classify the structural states based on the features. (no use of likelihood)

Feature Vector

\[ X = [ \omega_1, \omega_2, \omega_3, \omega_4, TMAC(H_0, \bar{\boldsymbol\phi}), MF(\bar{\boldsymbol\phi}, \bar{\boldsymbol\omega}) ] \]

where:

  • \(\omega_i\): is the \(i\)-th natural frequency.
  • \(TMAC(H_0, \bar{\boldsymbol\phi})\): is the Total Modal Assurance Criterion (TMAC) between the mode shapes of the undamaged state \(H_0\) and the observed mode shapes \(\bar{\boldsymbol\phi}\).
  • \(MF(\bar{\boldsymbol\phi}, \bar{\boldsymbol\omega})\): is the Modal Flexibility, which is a function of the observed mode shapes \(\bar{\boldsymbol\phi}\) and the observed natural frequencies \(\bar{\boldsymbol\omega}\).

Supervised learning - ML

General supervised learning model for classification: 1

\[ y = f(x, w) + \epsilon \]

  • \(y\): is the class (damage or no damage).
  • \(x\): is the feature vector.
  • \(w\): are the weights.
  • \(f(\cdot)\) is the function that maps the features to the class.
    • This is the decision rule “\(d(y)\)” that I have used erlier
  • \(\epsilon\): is the error term, which is there error in the model.

Training the models: 2

In general is the estimation of the weights \(w\) done by using Maximum Likelihood Estimation (MLE).

Note

Many problems do not have an closed-form solution for the MLE, and instead use optimization algorithms to find the weights that maximize the likelihood function.

ML Classification algorithms

Algorithms for classification:

Dataset

Data set for training the logistic regression model:

  • \(10^3\) samples for each sensor configuration
  • 50% class 0 (no damage)
  • 50% class 1 (damage)
    • equal number of each damage senario

Split:

  • 80% training data
  • 20% test data
# frequencyes and mode shapes
omega_all = load["natural_frequencies"]
phi_all = mode_shapes

# Coefficient of variation (CV)
CV_omega = 0.1
sigma_phi = 0.1

# Prior probabilities for the states (theta)
prior_theta = prior_theta

# Number of Samples 
N = 10**3

theta = OMA.theta
e_config_2 = e_configuration[e_configuration.sum(axis=1) >= 2]


theta = OMA.theta
bayes_risk_sensors = np.zeros(len(e_config_2))
cm = np.zeros((len(e_config_2), 2, 2), dtype=int) # confusion matrix for each sensor configuration

# Instantiate the BR class
A2 = BR(omega_all, phi_all, CV_omega, sigma_phi, prior_theta, N)


# improtance sampling prior for the states (theta) H_0 = 0.5 and H_1 = 0.5
prior_Importance = np.concatenate(([0.5], np.repeat(0.5 / (len(prior_theta) - 1), len(prior_theta) - 1)))

y_class = np.zeros(N*(len(e_config_2)))  # true class labels for each sensor configuration
X_features = np.zeros((N*(len(e_config_2)), 6))  # features: omega_bar, TMAC_H0, MF for each sensor configuration

# print(f"e_configuration[1:]:\n {e_configuration[1:]}")
for ei, S in enumerate(e_config_2):
    # Generate samples
    A2.MCS_modal(S, prior_Importance=prior_Importance)
    
    A2.omega_bar
    TMAC_H0 = np.array([A2.TMAC(A2.mask_and_normalize(A2.phi_H0, S), v) for v in A2.phi_bar])
    MF = np.array([A2.MF(phi, omega) for phi, omega in zip(A2.phi_bar, A2.omega_bar)])

    y_class_temp = np.array(A2.H_vec)  # true class labels (0 for H0, 1 for H1)
    X_features_temp = np.column_stack((A2.omega_bar, TMAC_H0, MF))

    y_class[ei*N:(ei+1)*N] = y_class_temp
    X_features[ei*N:(ei+1)*N, :] = X_features_temp
from sklearn.model_selection import train_test_split

# Split the dataset (80% train, 20% test)
X_train, X_test, y_train, y_test = train_test_split(
    X_features,
    y_class,
    test_size=0.2,      # 20% test
    random_state=42,    # ensures reproducibility
    shuffle=True        # ensures random split
)

print(f"Training set size: {len(X_train)} samples")
print(f"Test set size: {len(X_test)} samples\n")

print("Number of H0 and H1 samples in training set:", np.bincount(y_train.astype(int)))
print("Number of H0 and H1 samples in test set:", np.bincount(y_test.astype(int)))
Training set size: 8800 samples
Test set size: 2200 samples

Number of H0 and H1 samples in training set: [4269 4531]
Number of H0 and H1 samples in test set: [1114 1086]

PCA

Principal Component 1: [ 0.458  0.501  0.489  0.502  0.152 -0.16 ]
Principal Component 2: [-0.302  0.167  0.185  0.183  0.101  0.896]
Principal Component 3: [-0.083 -0.053 -0.101 -0.098  0.981 -0.088]
Principal Component 4: [ 0.801 -0.266 -0.364 -0.003  0.051  0.39 ]
Explained variance ratio: [0.41067787 0.17426818 0.16106897 0.09271319]
Total explained variance: 0.8387282067

\[ X_{\text{normalized}} = \frac{X - \mu}{\sigma} \]

Where:

  • \(X\): is the original feature vector.
  • \(\mu\): is the mean of the feature vector.
  • \(\sigma\): is the standard deviation of the feature vector.
Feature means: [ 0.95   6.14  17.319 31.201  0.986  1.769]
Feature stds: [0.18 0.93 2.62 4.6  0.01 4.39]

Logistic Regression

Observation vector:

\[ x = (x_1, ... , x_n) \]

Weighted sum: 1

\[ z = w_1 x_1 + ... + w_n x_n + b = \sum_{i=1}^{n} w_i x_i + b \]

Use logistic sigmoid function to turn it into a probability:

\[ \sigma(z)=\frac{1}{1+e^{-z}} \]

Thershold:

\[ \hat{y} = \begin{cases} 1 & \text{if } \sigma(z) \geq 0.5 \\ 0 & \text{if } \sigma(z) < 0.5 \end{cases} \]

from sklearn.linear_model import LogisticRegression

X, y = X_train, y_train

from sklearn import preprocessing
scaler = preprocessing.StandardScaler().fit(X)
print("Feature means before scaling:", np.round(scaler.mean_, 3))
print("Feature standard deviations before scaling:", np.round(scaler.scale_, 3))

X = scaler.transform(X)

# fit a logistic regression model to the data
clf = LogisticRegression(random_state=0).fit(X, y)

Feature means before scaling: [ 0.947  6.135 17.314 31.161  0.986  1.802]
Feature standard deviations before scaling: [0.183 0.918 2.614 4.621 0.014 4.885]


Classification Report:
               precision    recall  f1-score   support

         0.0       0.90      0.93      0.92      1114
         1.0       0.92      0.90      0.91      1086

    accuracy                           0.91      2200
   macro avg       0.91      0.91      0.91      2200
weighted avg       0.91      0.91      0.91      2200

Logistic Regression Coefficients: [[-0.987 -1.584 -1.601 -1.551 -0.639 -0.399]]
Logistic Regression Intercept: [0.785]

Analysis - Logistic regression

Note

  • this model performance is good for bad sensor configurations
  • But it is not as good for good sensor configurations, why?
    • I think that it is mostly based on the traning data.
from scipy.stats import multivariate_normal as mvn
from sklearn.metrics import confusion_matrix

# frequencyes and mode shapes
omega_all = load["natural_frequencies"]
phi_all = mode_shapes

# Coefficient of variation (CV)
CV_omega = 0.1
sigma_phi = 0.1

# Prior probabilities for the states (theta)
prior_theta = prior_theta

# Number of Samples 
N = 10**4

# must have an minimum of 2 sensors to be able to detect damage
e_config_2 = e_configuration[e_configuration.sum(axis=1) >= 2]


theta = OMA.theta
bayes_risk_sensors = np.zeros(len(e_config_2))
cm = np.zeros((len(e_config_2), 2, 2), dtype=int) # confusion matrix for each sensor configuration

# Instantiate the BR class
A2 = BR(omega_all, phi_all, CV_omega, sigma_phi, prior_theta, N)

# improtance sampling prior for the states (theta) H_0 = 0.5 and H_1 = 0.5
prior_H_Importance = np.array([0.5, 0.5])
prior_Importance = np.concatenate(([prior_H_Importance[0]], np.repeat(prior_H_Importance[1] / (len(prior_theta) - 1), len(prior_theta) - 1)))

print(f"Importance sampling prior probabilities for H: {prior_H_Importance}")
print(f"\nImportance sampling prior probabilities for θ: {np.round(prior_Importance, 4)}")

# print(f"e_configuration[1:]:\n {e_configuration[1:]}")
for ei, S in enumerate(e_config_2):
    # Generate samples
    A2.MCS_modal(S, prior_Importance=prior_Importance)

    A2.omega_bar
    TMAC_H0 = np.array([A2.TMAC(A2.mask_and_normalize(A2.phi_H0, S), v) for v in A2.phi_bar])
    MF = np.array([A2.MF(phi, omega) for phi, omega in zip(A2.phi_bar, A2.omega_bar)])

    # y_class_temp = np.array(A2.H_vec)  # true class labels (0 for H0, 1 for H1)
    X_features = np.column_stack((A2.omega_bar, TMAC_H0, MF))

    # d_vec = (Post1 > Post0).astype(int)
    X_test_scaled = (X_features  - scaler.mean_) / scaler.scale_
    d_vec = clf.predict(X_test_scaled)

    # Loss and Bayes Risk
    total_costs = prediction_cost[A2.H_vec, np.array(d_vec).astype(int)] + sensor_cost(len(e_config_2[sum(S)]))

    # Confusion matrix
    cm[ei] = confusion_matrix(A2.H_vec, d_vec)

    # Importance weights for the samples
    w = prior_H/prior_H_Importance
    importance_weights = w[A2.H_vec]  # Get the importance weights for each sample based on its true state H_vec


    # Estimation of the bayes risk using Monte Carlo.
    Bayes_Risk = np.mean(importance_weights * total_costs)

    print(f"Bayes Risk: {Bayes_Risk}")
    bayes_risk_sensors[ei] = Bayes_Risk
Importance sampling prior probabilities for H: [0.5 0.5]

Importance sampling prior probabilities for θ: [0.5    0.0333 0.0333 0.0333 0.0333 0.0333 0.0333 0.0333 0.0333 0.0333
 0.0333 0.0333 0.0333 0.0333 0.0333 0.0333]
Bayes Risk: 29.130507999999995
Bayes Risk: 27.343752000000002
Bayes Risk: 22.956972
Bayes Risk: 20.93484
Bayes Risk: 25.122072
Bayes Risk: 25.169428
Bayes Risk: 21.5434
Bayes Risk: 28.200376000000002
Bayes Risk: 20.380892000000003
Bayes Risk: 19.34818
Bayes Risk: 18.221691999999997

Minimum Bayes Risk: 18.221691999999997
Best sensor configuration: [1 1 1 1] (index: 10)

Questions ?

Error In OMA?

What is ralistic values of variaction in the error in OMA?

  • frequensy
  • mode shapes

Note

  • \(\phi\) is assumed to be uniform, between -1 and 1
  • \(\phi\) is normilized, so it has at least one value of 1 or -1.
  • \(\phi\) is there for in this cass somewhat dependent on the lengh of the mode shape vector.

Mean values of |phi|: 0.498689392669938
Number of samples: 1000

Mean values of |phi|: 0.5073268141520177
Number of samples: 1000
Mode shapes 0:
[[ 0.093  0.328  0.647  1.   ]
 [ 0.505  1.     0.544 -0.726]
 [ 1.     0.335 -0.972  0.427]
 [ 1.    -0.968  0.617 -0.174]]

Sum of the absolute values of the first mode shape across all samples: 
[2.59767468 2.63070467 2.78064049 2.32723371]

Mean of the absolute values of the first mode shape across all samples: 
[0.64941867 0.65767617 0.69516012 0.58180843]

damage and material properties

Question ?

  • How does damage affect the material properties of the structure?
  • How does damage affect the probability of structural failure?

:::

damage level?

Question ?

What damage level do my model need to be able to detect?

Idear 1

  • Introduce probaility of failure function
  • Introduce probability of stiffness change function

Probability of failure given the damage level and the structural state:

\[ P(\text{failure} \mid \Delta E, \theta) , \text{ continuous function} \]

Probability of stiffness change given structural states: \[ P(\Delta E \mid H) , \text{ continuous function} \]

What will we use them for?

  • Update the cost matrix, so that the cost of a false negative is given as the probability of failure times the expected cost of failure

Assumptions / statements

  • The estimation of thes probalistic function requires an separate simulation study.
    • make assumptions about the functions ?

Reformulation of model to continuous damage levels:

  • Gives no significant changes in the logistic regression model, only in how the data is generated?
  • change the formulation of the likelihood function. (Need to think)

Whats next?

My next work suggestions

  • Week +1: Work on model capability.
  • Week +2: Bigger FEM model + more realistic Priors and cost functions.

.

  • 3 structural casse studies
    • Hospital building (hige cost of FN)
    • ~ (hige cost of FP)
    • Low cost of both FN but high relative cost of FP
  • bigge FEM model