Case Study

Cantilever beam

Thyge Vinther Ludvigsen, s203591

15 Apr, 2026

Introduction

Goal of the case study

What I test in this case study:

simple link function/ Error function.
overall workflow.
Evaluation methods of the Bayesian risk

Note:

compute the full sensor solution space
- no optimization, no sensor selection.

Methodology

Chosen Features

Modal properties:

2 first natural frequencies

Features:

2 first natural frequencies

\[ 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.

Prior distribution

Theory
Venn diagram
Code

No Damage and Damage: \[ P(H_0) = P(\theta_0) = 0.90 \]

\[ P(H_1) = 1 - P(H_0) \]

Assume that all damage states with damage are equally likely

\[ P(\theta_i \mid H_1) = \frac{1}{N} \quad \text{for } i = 1, 2, ..., K \]

This gives. (Remember that \(P(\theta_i) = P(\theta_i, H_1)\) since \(\theta_i\) can only occur if \(H_1\) is true) \[ P(\theta_i) = P(\theta_i, H_1) = P(H_1) P(\theta_i \mid H_1) = P(H_1) \frac{1}{N} \quad \text{for } i = 1, 2, ..., K \]

where \(N\) is the number of damage states considered.

\[ N = \sum \frac{n!}{k!(n-k)!} = 2^n = 2^4 = 16 \]

n is the number of elements
k is the number of damaged elements
the binomial coefficient \(\frac{n!}{k!(n-k)!}\) gives the number of ways to choose k damaged elements from n total elements.

\(P(\theta_i) = P(\theta_i, H_1)\) for \(i = 1, 2, ..., K\) because \(\theta_i\) can only occur if \(H_1\) is true.

P_H0 = 0.9
P_H1 = 1-P_H0

prior = np.zeros(16)
prior[0] = P_H0
prior[1:] = P_H1/len(prior[1:])

print(f"prior: {prior}")

prior: [0.9        0.00666667 0.00666667 0.00666667 0.00666667 0.00666667
 0.00666667 0.00666667 0.00666667 0.00666667 0.00666667 0.00666667
 0.00666667 0.00666667 0.00666667 0.00666667]

Error function \(\epsilon(e)\)

Note: this error function is only for the estimation of the natural frequencies.

\[ \epsilon_{\omega}(e) = \epsilon_0(e) + ... + \epsilon_n(e) \]

Each error term explains a different source of error in the estimation of the natural frequencies. I assume that the function above can be constructed in a way where all the error terms are independent, and that the total error can be calculated as the sum of the individual error terms. (This is not True, example of signal to noise ratio and the effect of temperature to the stiffness of the structure.)

Error form Sensor signal noise

Model
\(\bar{X} \sim \mathcal{N}\left(\mu, \frac{1}{n^2}\sum_{i=1}^n \sigma_i^2\right)\)

Assumptions

I assume that the natural frequencies is estimated as the mean of all the OMA estimations of the natural frequencies for all the sensor locations.
- \(\hat{\omega} = \frac{1}{n} \sum_{i=1}^n \omega_i\) where \(\omega_i\) is the estimation of the natural frequencies from sensor location \(i\).
I assume that the error from environmental effects to be independent of the sensor noise.
I assume that the error from the sensor noise in the prediction of the natural frequencies is normally distributed, and that the variance of this error is a function of the sensor locations and the sensor noise level.

This assumptions leads to.

\[ \epsilon_{vibration} \sim \mathcal{N}(0, \frac{1}{\# e_{vibration}^2}\sum_{i=1}^{\# e_{vibration}} \sigma(e_{vibration, i})^2) \]

Where:

\(\epsilon_{vibration}\) is the error from sensor noise in the estimation of the natural frequencies.
\(\#\) is the count operator.
\(\sigma(e_{vibration, i})\) is the standard deviation of the error for an given sensor \(i\).

Important

We can have an increasing error with more sensors, if the error from the sensor noise is high.

We have \(n\) independent normal distributions where we make 1 draw from each distribution.

\(X_1, X_2, \dots, X_n \text{ independent, } X_i \sim \mathcal{N}(\mu, \sigma_i^2)\)

we compute the sample mean.

\(\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i\)

now do we find the distribution of \(\bar{X}\), it is know that the sample mean of a normal distribution is also normally distributed.

\[ \mathbb{E}[\bar{X}] = \mathbb{E}\left[\frac{1}{n}\sum_{i=1}^n X_i\right] = \frac{1}{n}\sum_{i=1}^n \mathbb{E}[X_i] = \frac{1}{n}(n\mu) = \mu \]

the draws are independent so the variance of the sample mean is:

\[ \mathrm{Var}(\bar{X}) = \mathrm{Var}\left(\frac{1}{n}\sum_{i=1}^n X_i\right) = \frac{1}{n^2} \sum_{i=1}^n \sigma_i^2 \]

note: \(1/n^2\) is from the scaling rule of the variance \(Var(aX + b) = a^2 Var(X)\), it can easily be seen by moving \(1/n\) inside the sum.

final result:

\[ \bar{X} \sim \mathcal{N}\left(\mu, \frac{1}{n^2}\sum_{i=1}^n \sigma_i^2\right) \]

Linear combination of non independent random variables

Std, \(\sigma(e_{vibration, i})\)

Simplest model
Variables
More complex model ??

Minimum require

The model need to weight the positions of the sensors.
- use mode shapes
and the number of sensors.
- this is done in the variance of the distribution function for the mean error \(\frac{1}{n^2}\sum_{i=1}^n \sigma_i^2\)

\[ \sigma(e_{vibration, i}) = \omega_i \cdot 5\% \cdot f(\phi_i(e_{vibration, i})) \]

where:

\(\omega_i\) is the natural frequency of mode \(i\).
\(5\%\) is the assumed sensor coefficient of variation (CV)
\(f(\phi_i(e_{vibration, i})) = 10 - 9 \cdot Abs(\phi_i(e_{vibration, i}))\) is a function that weights the sensor location to noise experienced in the signal.
- \(\phi_i(e_{vibration, i})\) is the value of the normalized mode shape at the sensor location \(e_{vibration, i}\) for mode \(i\). this is a value between 0 and 1, where 0 is no vibration and 1 is the maximum vibration.

Modal variables:

Natural frequencies: \(\omega_1, \omega_2\)
Mode shapes: \(\phi_1, \phi_2\)
Participation factors: \(\Gamma_1, \Gamma_2\)
- Needs to be the participation in the response of tee structure for an given ambient excitation.
- Assume it to be the same as mass participation factors ?

Model setup variables:

Sensor configuration: \(e_{dva}\)
length of the signal, \(T\)
- number of cycles captured in the signal, \(n_c\)
- number of samples captured in 1 cycle for an given mode, \(n_s\)
- sampling frequency of the signal, \(f_s\)

Ambient excitation variables:

Amplitude of the Mode shape given the ambient excitation. \(A_{i}\)
The Ambient excitation effect the participation factors for each mode, it is dependent on the frequency content of the ambient excitation.

Sensor variables:

Signal to noise ratio (SNR) of the sensor
Sampling frequency of the sensor, Hz \(f_s\)
Sensor type (accelerometer, velocity sensor, displacement sensor)
Sensor location (\(e_{dva}\))

A more complex model need to both include almost all of the Variables. Time and testing is need to formulate this model

Code exampledef sigma_vibration(omega, phi, node_number):
    # assumed sensor coefficient of variation (CV)
    CV = 0.05
    # calculate vibration standard deviation
    sigma = omega * CV * (10 - 9 * abs(phi[node_number]))
    return sigma

print(f"omega_1 = {omega_1:.2f} s^-1, mode shape phi_1 = {phi_1}")
print(f"omega_2 = {omega_2:.2f} s^-1, mode shape phi_2 = {phi_2}")

node_number = 4
print(f"node {node_number}: sigma_vibration for mode 1 = {sigma_vibration(omega_1, phi_1, node_number):.4f}")
print(f"node {node_number}: sigma_vibration for mode 2 = {sigma_vibration(omega_2, phi_2, node_number):.4f}")

omega_1 = 1.04 s^-1, mode shape phi_1 = [0.         0.09252781 0.32807894 0.64695444 1.        ]
omega_2 = 6.72 s^-1, mode shape phi_2 = [ 0.          0.50514687  1.          0.54397217 -0.72646916]
node 4: sigma_vibration for mode 1 = 0.0522
node 4: sigma_vibration for mode 2 = 1.1635

Setup 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 =  list(itertools.product([0, 1], repeat=len(e_sensor))) # all combinations of 0 and 1 for each sensor node
e_configuration = e_configuration* e_sensor 
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 4]
Configuration 2 given in nodes with sensors: [0 0 3 0]
Configuration 3 given in nodes with sensors: [0 0 3 4]
Configuration 4 given in nodes with sensors: [0 2 0 0]
Configuration 5 given in nodes with sensors: [0 2 0 4]
Configuration 6 given in nodes with sensors: [0 2 3 0]
Configuration 7 given in nodes with sensors: [0 2 3 4]
Configuration 8 given in nodes with sensors: [1 0 0 0]
Configuration 9 given in nodes with sensors: [1 0 0 4]
Configuration 10 given in nodes with sensors: [1 0 3 0]
Configuration 11 given in nodes with sensors: [1 0 3 4]
Configuration 12 given in nodes with sensors: [1 2 0 0]
Configuration 13 given in nodes with sensors: [1 2 0 4]
Configuration 14 given in nodes with sensors: [1 2 3 0]
Configuration 15 given in nodes with sensors: [1 2 3 4]

Errors for all sensor configurations

\[ f(y|\theta) \sim \mathcal{N}(\mu, \Sigma^2), \text{ where } \Sigma^2 = \frac{1}{n^2}\sum_{i=1}^n \sigma_i^2 \]

2 first natural frequencies

Undamaged case:
Natural frequencies (Hz = s^-1): [1.04366091 6.72223044]

Damaged cases:
Natural frequencies (Hz = s^-1):
[[1.0126304  4.63016511]
 [0.80548911 3.69172521]
 [0.78872241 3.37286305]
 [0.56004051 4.87740021]
 [0.55489562 4.11077885]
 [0.5120402  3.3640654 ]
 [0.50788611 3.08229423]
 [0.39753106 4.75750033]
 [0.39595636 3.78110917]
 [0.38094819 2.57381063]
 [0.37948544 2.37804179]
 [0.34212249 2.92025964]
 [0.34106067 2.71671264]
 [0.33103083 2.23314855]
 [0.33003456 2.12575592]]

Analysis: Monete Carlo

Monte Carlo simulation (MCS)

Workflow
Big-O

Sudocode for the MCS:

… Othere code

for any sensor configuration \(e_i\):
N simulations:
- Sample a damage state \(\theta_i\) from the prior distribution \(P(\theta)\)
- Simulate the sensor measurements \(y_i\) given the damage state \(\theta_i\) and the sensor configuration \(e_i\)
- compute the pdf of \(y\) given \(\theta_i\) and \(e_i\), \(f(y|\theta_i)\)
- Compute the posterior distribution \(P(H|y_i)\) using Bayes’ theorem
- Compute the expected cost for the sensor configuration \(e_i\) using the posterior distribution and the cost function

… Othere code

Symbol	Values	Meaning
\[N\]	\(\sim 10^{3-7}\)	Number of outer Monte Carlo samples (Bayes risk estimation)
\[H\]	\(2\)	Number of structural hypotheses \[H_i\]
\[\Theta\]	\(2^{(\text{number of elements in the structure})}\)	Number of \[\theta\]-samples used to approximate \[f(y\mid H)\]
\[G\]	small	Cost of evaluating the forward model \[g(\theta,e)\] * (\(g(x(\theta,e))\) is expensive)
\[Y\]	2-100	Dimension of feature vector \[y\]

\[ O(N \cdot [ G + H \cdot \Theta \cdot Y^{(*)}]) \]

where:

The term \(Y^{(*)}\) is the cost of generating the error for \(y\).
- if all feature error are independent the is \(Y^{(*)} = Y\) since we can generate the error for each feature independently.
- if the errors for the features are dependent but \(\Sigma_{i,j} = \Sigma_{j,i}\) then \(Y^{(*)} = Y^2\).
- if the errors for the features are dependent and \(\Sigma_{i,j} \neq \Sigma_{j,i}\) then \(Y^{(*)} = Y^3\).

How to reduce the cost of the MCS:

\[ O(N \cdot [ G + H \cdot \Theta_{\text{surrogate}} \cdot Y^{1}] + \text{Cost of surrogate model}) \]

The minimum cost is by keeping the feature error independent
Make an surrogate model of \(f(y\mid H)\) outside the MCS simulations, use this surrogate model to approximate \(f(y\mid H)\) within the MCS simulations.

Illustration - MCS

Level 1 SHM system

flowchart TB
    subgraph SensorDesign["Sensor configurations: e    ."]
        direction TB
        style SensorDesign fill:#e8f5e9,stroke:#2e7d32,stroke-width:3px,color:#1b5e20
        
        subgraph DamageNoDamage["Structural categories: H    ."]
            direction TB
            style DamageNoDamage fill:#fff9c4,stroke:#f57f17,stroke-width:3px,color:#e65100
            
            subgraph Measurement["Feature space: y    ."]
                direction TB
                style Measurement fill:#fff3e0,stroke:#e65100,stroke-width:3px,color:#bf360c

                A["Features<br/>y = g(θ, e) + ε(e)"]
                A --> B

                subgraph StructuralState["Structural states: θ    ."]
                    direction TB
                    style StructuralState fill:#e3f2fd,stroke:#1565c0,stroke-width:3px,color:#0d47a1
                
                    B["Likelihood<br/>f(y | H) = ∑<sub>θ</sub> f(y | θ) P(θ | H)"]
                    end

            C["Posterior<br/>P(H | y)"]
            D["Decision Rule<br/>d = max P(H | y)"]
            E["Frequentist Risk<br/> R(d) = ∫ f(y | H) L(H, d) dy"]

            B --> C
            C --> D
            D --> E
            end

            F["Bayes Risk<br/> Ψ(e) = ∑<sub>H</sub>  P(H) R(d)"]
            E --> F
        end

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

Monte Carlo Simulation (MCS)

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


            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

MCS - Code

Question / Note

\(f(y | \theta)\) is not an real probability density, it is just the pdf at the given observation \(y\) for a given state \(\theta\). the probability of an continuous random variable taking a specific value is 0.

"""Utilize vectorization and reduce the calls to np.random.choice and np.random.multivariate_normal."""
import numpy as np
from scipy.stats import multivariate_normal
np.random.seed(42)

# Number of Monte Carlo samples
N = 10**4

# True state features, damage states, initialization of bayes risk array
y_start = load["natural_frequencies"][0:, :2]
theta = OMA.theta
bayes_risk_sensors = np.zeros(len(e_configuration))

# Loop over sensor configurations
for ei in range(len(e_configuration)):
    # sensor configuration + covariance matrix
    e = e_configuration[ei]
    cov = np.diag(sigma[:, ei]**2)
    
    # Generate all states and H labels at once: theta ~ P(theta)
    states = np.random.choice(range(len(prior)), size=N, p=prior)
    H_vec = (states > 0).astype(int)  # 0 if state == 0, else 1
    
    # Generate all errors and observations (y) at once: y = g(theta) + epsilon, epsilon ~ P(epsilon|e, theta)
    epsilon_vec = np.random.multivariate_normal(mean=np.zeros(2), cov=cov, size=N)
    y_vec = np.maximum(y_start[states] + epsilon_vec, 0)
    
    # Vectorized Likelihood calculation:  pdf = f(y|theta_i)
    pdf_epsilon_all = np.zeros((N, len(theta)))
    for i in range(len(theta)):
        pdf_epsilon_all[:, i] = multivariate_normal.pdf(y_vec - y_start[i], mean=np.zeros(2), cov=cov)
        
    # Likelihood: f(y|H) = sum_i f(y|theta_i) * P(theta_i|H)
    Likelihood_H0 = pdf_epsilon_all[:, 0] * prior[0] / P_H0 # 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[1:] / P_H1), axis=1)
    
    # Posterior probabilities: P(H|y) ∝ f(y|H) * P(H)
    Post0 = Likelihood_H0 * P_H0
    Post1 = Likelihood_H1 * P_H1

    # Decisions rule: maximimum a posteriori (MAP) 
    d_vec = (Post1 > Post0).astype(int)
    
    # Loss and Bayes Risk
    total_costs = prediction_cost[H_vec, d_vec] + sensor_cost(len(e_configuration[ei]))
    
    # 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: 66.115
Bayes Risk: 298.9
Bayes Risk: 110.58
Bayes Risk: 18.955
Bayes Risk: 61.105
Bayes Risk: 129.975
Bayes Risk: 82.1
Bayes Risk: 520.36
Bayes Risk: 195.61
Bayes Risk: 249.87
Bayes Risk: 128.57
Bayes Risk: 152.575
Bayes Risk: 79.01
Bayes Risk: 108.23
Bayes Risk: 65.59

Results - MCS

Risk for Uninformed system

Assuming that.

the decision is always \(d_0\) (no damage predicted).
we can use the prior distribution to calculate the risk for this decision.

\[ R(d_0) = C_\text{prediction}(H_0, d_0) \cdot P(H_0) + C_\text{prediction}(H_1, d_0) \cdot P(H_1) \]

# C [H, d] = the prediction cost
prediction_cost[0, 0] * P_H0 + prediction_cost[1, 0] * P_H1

np.float64(999.9999999999998)

Risk for informed system

Minimum Bayes Risk: 18.955
Best sensor configuration: (np.int64(2),) (index: 3)

The results here are highly dependent on:

the error function \(\epsilon(e)\)
the cost function \(C(H, d)\)
the prior distribution \(P(H)\) and the distribution of the damage states \(P(\theta)\)
the chosen degree of damage in the damage states \(\theta\).

Analysis: Mean value

Use the mean value approximation to evaluate the Bayes risk.

Mean value approximation / Point estimate method

Idea:

Simplify all stochastic variables to their mean value, and then evaluate the Bayes risk for this simplified system.

Down sides:

it is not accurate, and ignores most of the uncertainty in the system.

Code example"""Mean value approximation"""
import numpy as np
from scipy.stats import multivariate_normal

# True state features, damage states, initialization of bayes risk array
y_start = load["natural_frequencies"][0:, :2]
theta = OMA.theta
bayes_risk_sensors = np.zeros(len(e_configuration))

# Loop over sensor configurations
for ei in range(len(e_configuration)):
    # sensor configuration + covariance matrix
    e = e_configuration[ei]
    cov = np.diag(sigma[:, ei]**2)
    
    # Loop over states (damage cases): theta
    for theta_i in range(len(theta)):
        # print(f"theta_i: {theta_i}, damage case: {theta[theta_i]}")
        H_vec = int(theta_i > 0)  # 0 if state == 0, else 1
    
        # makes y whit 0 error
        y_vec = np.maximum(y_start[theta_i], 0)
        
        # Vectorized Likelihood calculation:  pdf = f(y|theta_i)
        pdf_epsilon_all = np.zeros( len(theta))
        for i in range(len(theta)):
            pdf_epsilon_all[i] = multivariate_normal.pdf(y_vec - y_start[i], mean=np.zeros(2), cov=cov)
            
        # Likelihood: f(y|H) = sum_i f(y|theta_i) * P(theta_i|H)
        Likelihood_H0 = pdf_epsilon_all[0] * prior[0] / P_H0 # 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[1:] / P_H1), axis=0)
        
        # Posterior probabilities: P(H|y) ∝ f(y|H) * P(H)
        Post0 = Likelihood_H0 * P_H0
        Post1 = Likelihood_H1 * P_H1

        # Decisions rule: maximimum a posteriori (MAP) 
        d_vec = (Post1 > Post0).astype(int)
        
        # Loss and Bayes Risk
        total_costs = prediction_cost[H_vec, d_vec] + sensor_cost(len(e_configuration[ei]))
        
        # 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: 106.0
Bayes Risk: 106.0
Bayes Risk: 107.0
Bayes Risk: 106.0
Bayes Risk: 107.0
Bayes Risk: 107.0
Bayes Risk: 108.0
Bayes Risk: 106.0
Bayes Risk: 107.0
Bayes Risk: 107.0
Bayes Risk: 108.0
Bayes Risk: 107.0
Bayes Risk: 108.0
Bayes Risk: 108.0
Bayes Risk: 109.0

Results - Mean value approximation

Minimum Bayes Risk: 106.0
Best sensor configuration: (np.int64(4),) (index: 0)

Analysis: Numerical integration

Pappers:

Method is proposed in the paper:

“A univariate dimension-reduction method for multi-dimensional integration in stochastic mechanics” (Rahman and Xu, 2004, p. 393) (pdf)

The method is used in the paper:

“An optimal sensor placement design framework for structural health monitoring using Bayes risk” (Yang et al., 2022, p. 1) (pdf)

Cost of Evaluating the Bayes risk

Evaluate the 3 methods for evaluating the Bayes risk:

Time comparison of the 3 methods

Monte Carlo simulation, whit \(N=10^4\): 0.6 seconds
Mean value approximation: 0.1 seconds
Numerical integration:

Sensativity analysis

Feacture states

Minimum Bayes Risk: 951.0
Best sensor configuration: (np.int64(3), np.int64(4)) (index: 2)

Minimum Bayes Risk: 882.13
Best sensor configuration: (np.int64(4),) (index: 0)

Minimum Bayes Risk: 601.64
Best sensor configuration: (np.int64(4),) (index: 0)

Minimum Bayes Risk: 439.83
Best sensor configuration: (np.int64(4),) (index: 0)

Minimum Bayes Risk: 301.73
Best sensor configuration: (np.int64(2),) (index: 3)

Minimum Bayes Risk: 224.005
Best sensor configuration: (np.int64(2),) (index: 3)

Minimum Bayes Risk: 110.02
Best sensor configuration: (np.int64(2),) (index: 3)

Minimum Bayes Risk: 65.17
Best sensor configuration: (np.int64(2),) (index: 3)

Minimum Bayes Risk: 20.01
Best sensor configuration: (np.int64(2),) (index: 3)

Cost function

FT = 0

cost = np.array([[0,   0],  # TN, FP
                 [10000, 100]   # FN, TP
                 ])  
sensativity_cost_function(cost)

FT = 1000

TP = 0

cost = np.array([[0,   50],  # TN, FP
                 [10000, 0]   # FN, TP
                 ])  
sensativity_cost_function(cost)

TP = 10000

FN = 100

cost = np.array([[0,   50],  # TN, FP
                 [100, 100]   # FN, TP
                 ])  
sensativity_cost_function(cost)

FN = 100000

Prior probability distribution

Cost function (prediction_cost): 
[[    0    50]
 [10000   100]]

Minimum Bayes Risk: 18.955
Best sensor configuration: (np.int64(2),) (index: 3)

Minimum Bayes Risk: 7.18
Best sensor configuration: (np.int64(2),) (index: 3)

Minimum Bayes Risk: 6.1
Best sensor configuration: (np.int64(2),) (index: 3)

Minimum Bayes Risk: 6.0
Best sensor configuration: (np.int64(3),) (index: 1)

Minimum Bayes Risk: 6.0
Best sensor configuration: (np.int64(3),) (index: 1)

More realistic case study + number of simulations

prediction_cost = np.array([[0,   100_000],  # TN, FP
                 [10_000_000, 1_000_000]   # FN, TP
                 ])

P_H1 = 10**(-4)
P_H0 = 1 - P_H1
prior = np.array([P_H0] + [P_H1/(len(theta)-1)]*(len(theta)-1))

Minimum Bayes Risk: 6.0
Best sensor configuration: (np.int64(4),) (index: 0)

Minimum Bayes Risk: 6.0
Best sensor configuration: (np.int64(2),) (index: 3)

Minimum Bayes Risk: 66.0
Best sensor configuration: (np.int64(2),) (index: 3)

Minimum Bayes Risk: 144.1
Best sensor configuration: (np.int64(2),) (index: 3)

What’s next

Discussion

Main direction
Resmblance to the real world
Capabilities
Case study

The model’s resemblance to the real world (or the assumptions underlying it).
The model’s capabilities based on those assumptions.
Case study results with realistic data.

Prior probability distribution
Cost function
Error function link to sensor configuration (major assumption!)
- Make an better link function.
Include more stochasticity in the model: Temperature, Mass changes, unknown base model,…
Continuous damage states

Feature selection methods (e.g., PCA)
Inclusion of additional features (e.g., MAC)
Optimization algorithms for sensor placement
Short-term memory: Incorporate the last 10 observations instead of relying on a single observation
Improved evaluation of Bayes risk:
- Importance sampling for Monte Carlo Simulation (MCS)
- Numerical integration (e.g., method by Rahman and Xu, 2004)
- Surrogate models (???)

Resarch realistic data for.
- cost function
- prior distribution
- error function / link function (?)

My Resening

My Resening
Resmblance to the real world
Capabilities
Case study

I have chosen to split the work for next time into 4 main class:

Do it now
Can be don leater
Can be chosen as an Asumption or Simplification
Other…

Can be don leter:

Prior probability distribution
Cost function

Asumption / Simplification:

Include more stochasticity in the model: Temperature, Mass changes, unknown base model,…
Continuous damage states

Other:

Error function link to sensor configuration (major assumption!)
- Make an better link function.

Do it now:

Feature selection methods (e.g., PCA)
Inclusion of additional features (e.g., MAC)
Optimization algorithms for sensor placement
Improved evaluation of Bayes risk:
- Importance sampling for Monte Carlo Simulation (MCS)
- Numerical integration (e.g., method by Rahman and Xu, 2004)
- Surrogate models (???)

Can be done later:

Short-term memory: Incorporate the last 10 observations instead of relying on a single observation

Can be done later:

Research realistic data for.
- cost function
- prior distribution
- error function / link function (?)
Make assumptions for Geometry and material properties.

Time needed for tasks

Improved evaluation of Bayes risk: (1-3 weeks)

Importance sampling (MCS): 1-3 days
Numerical integration: 1-2 weeks (seams to be quite complex to understand)
Surrogate models (???): ???

Inclusion of additional features: ( ???? )

Formulation of new error function
- What assumptions to use?

Optimization algorithms for sensor placement: (2-7 days)

Make an new structure whit more sensor configurations: 1-3 days
- I already have some papers there refere to some algorithms.
Research, test and implement optimization algorithms: 1-3 days

Feature selection methods (e.g., PCA): (1 week)

Research literature on feature selection methods: +3 days
Implement and test feature selection methods: 1-3 days

Questions

How to best link the error/stochasticity of the sensor configuration to features???

Approache 1:

Make Asumption for all the tings???

Approache 2:

Create data, and make assumptions based on this data.
- Setup model stochasticity
- Simulate sensor data
- OMA on sensor data
- Calculate features
- Make assumption based on the link between the sensor configuration and the error in the features.
note: this problem is not straight forward and have many potentially time consuming steps.

Approache 3:

Research literature: I have looked but i can look again.