i want to plot fragility curve stress vrs probability of failure and mean speed vrs probability of failure no of data is 48500

% Load the data
filename = 'vehicle_data7.xlsx';  % Data file
data = readtable(filename);       % Load the data into a table

% Beam properties
L = 40;                    % Beam length (m)
h = 2;                     % Beam height (m)
allowable_stress = 240e6;   % Allowable stress in Pa (240 MPa)
g = 9.81;                  % Gravitational acceleration (m/s^2)

% Constants
rho_steel = 7850;           % Density of steel (kg/m^3)

% Number of simulations (equal to number of rows in dataset)
num_simulations = height(data);

% Initialize counters for failures
failure_count = 0;

% Monte Carlo simulation loop
for sim = 1:num_simulations
    % Extract vehicle data
    vehicle_type = data.VehicleType{sim}; % Vehicle type (e.g., sedan, bus, truck)
    weight = data.Weight_kN(sim) * 1e3;   % Convert from kN to N
    speed = data.Speed_km_h(sim) / 3.6;   % Convert from km/h to m/s
    spacing = data.Spacing_m(sim);        % Spacing (m)
    E_veh = data.Modulus_of_Elasticity_GPa(sim) * 1e9; % Convert from GPa to Pa
    I_veh = data.Moment_of_Inertia_m4(sim); % Moment of inertia (m^4)
    damping_ratio = data.Damping_Ratio(sim); % Damping ratio from dataset

    % Beam's cross-sectional area and mass per unit length
    b = I_veh * 12 / h^3;                 % Calculate beam width from moment of inertia
    A = b * h;                            % Cross-sectional area (m^2)
    m = rho_steel * A;                    % Mass per unit length (kg/m)

    % Modal mass (corrected formula)
    m_n = (m * L) / 2;                    % Modal mass (kg)

    % Natural frequency for the first mode (omega_1)
    omega_1 = sqrt((E_veh * I_veh) / (m_n * L^4));  % Natural frequency for first mode (rad/s)

    % Stiffness for the first mode
    k_1 = m_n * omega_1^2;               % Modal stiffness (N/m)

    % Critical velocity (m/s)
    v_c = omega_1 * (L / pi);            % Critical velocity

    % Dynamic Amplification Factor (DAF) calculation
    DAF = 1 / sqrt((1 - (speed / v_c)^2)^2 + (2 * damping_ratio * (speed / v_c))^2);


    % Mode shape function for the first mode (phi_1)
    phi_1 = @(x) sin(pi * x / L);        % Mode shape for first mode (phi_1)
    F_1 = weight * phi_1(L / 2);         % Modal force due to moving load (N)

    % Maximum displacement
    w_max = (F_1 / k_1) * DAF;          % Maximum dynamic displacement (m)

    % Maximum stress (sigma_max)
    phi_1_dd = @(x) -(pi^2 / L^2) * sin(pi * x / L); % Second derivative of phi_1
    z = h / 2;                           % Distance from neutral axis to outer fiber (m)
    M_max = -E_veh * I_veh * phi_1_dd(L / 2) * w_max; % Maximum bending moment (N.m)
    sigma_max = M_max * z / I_veh;       % Maximum dynamic stress (Pa)

    % Compare sigma_max with allowable stress to determine failure probability
    if sigma_max > allowable_stress
        failure_count = failure_count + 1;  % Increment failure counter
    end
end

% Calculate failure probability
failure_probability = failure_count / num_simulations;

% Display results
fprintf('Number of Simulations: %d\n', num_simulations);
fprintf('Number of Failures: %d\n', failure_count);
fprintf('Failure Probability: %.4f\n', failure_probability);

Answer 1

[Walter Roberson]

Your failure probability is a scalar. mean speed could be calculated as mean(data.Speed_km_h) but that would be a scalar. Plotting a scalar vs a scalar is a completely uninteresting plot.

You could conceivably

    if sigma_max > allowable_stress
        failure_count = failure_count + 1;  % Increment failure counter
        isfailure(sim) = true;
    else
        isfailure(sim) = false;
    end

and then plot

    scatter(data.Speed_km_h, isfailure)

I have to wonder, though, whether you should be doing something like discretize() Speed_km_h and accumarray() isfailure against discretized speed, and then bar plot the result.

[AL]

yes but how can ???actually i want to plot fragility curve but i cant..and i also confused whether i plot probability of failure vrs stress or probability of failure vrs mean velocity ??please suggest me how can i solve ??

[Walter Roberson]

It is not clear to me how fragility is to be calculated.

Probability of failure is a scalar. You cannot plot a scalar against a vector... not unless you are expecting a horizontal or vertical line as the only result of the plot.

The one way you could plot probability of failure meaningfully is if you had a number of different files to read in, and you wished to calculate probability of failure of each of them.

I think you would do much better to plot a histogram of failure count relative to speed.

[AL]

this is my random data generation code .actually i want to plot fragility curve probability of failure vrs mean velocity could you give me suggestion how can i developed my code

clc;
clear all;
close all;
% Define vehicle types
vehicle_types = {'Sedan', 'Bus', 'Truck'};
% Parameters: mean and standard deviation for weight, speed, and spacing
weight_params = [12, 2; 150, 20; 100, 50];  % Rows Sedan, Bus, Truck
speed_params = [50, 20; 50, 20; 50, 20];  % Rows: Sedan, Bus, Truck
spacing_params = [50, 20; 50, 20; 50, 20];  % Rows: Sedan, Bus, Truck
% New parameters: modulus of elasticity, moment of inertia, damping ratio
elasticity_params = [200, 5; 200, 5; 200, 5];  % Mean and std for all
inertia_params = [0.02, 0.001; 0.02, 0.001; 0.02, 0.001];  % Mean and std for all
damping_params = [0.01, 0.0005; 0.01, 0.0005; 0.01, 0.0005];  % Mean and std for all
% Number of examples for each vehicle type
num_examples = [40000, 500, 8000];  % Sedan, Bus, Truck
% Initialize storage for the dataset
random_data = [];
% Generate random data for each vehicle type
for i = 1:length(vehicle_types)
    % Extract mean and standard deviation for weight
    mean_weight = weight_params(i, 1);
    std_weight = weight_params(i, 2);
    % Check if the vehicle type is Truck to generate log-normal weights
    if strcmp(vehicle_types{i}, 'Truck')
        % Compute log-normal parameters
        mu = log((mean_weight^2) / sqrt(std_weight^2 + mean_weight^2));
        sigma = sqrt(log(1 + (std_weight^2) / (mean_weight^2)));
        % Generate log-normal distributed weights for Trucks
        weights = lognrnd(mu, sigma, [num_examples(i), 1]);
    else
        % Generate normal distributed weights for Sedan and Bus
        weights = normrnd(mean_weight, std_weight, [num_examples(i), 1]);
    end
    % Generate normal distributed speed and spacing for all vehicle types
    speeds = normrnd(speed_params(i, 1), speed_params(i, 2), [num_examples(i), 1]);
    spacings = normrnd(spacing_params(i, 1), spacing_params(i, 2), [num_examples(i), 1]);
    % Generate new parameters
    elasticity = normrnd(elasticity_params(i, 1), elasticity_params(i, 2), [num_examples(i), 1]);
    inertia = normrnd(inertia_params(i, 1), inertia_params(i, 2), [num_examples(i), 1]);
    damping = normrnd(damping_params(i, 1), damping_params(i, 2), [num_examples(i), 1]);
    % Create a table for this vehicle type
    vehicle_data = table(repmat(vehicle_types(i), num_examples(i), 1), ...
        weights, speeds, spacings, elasticity, inertia, damping, ...
        'VariableNames', {'VehicleType', 'Weight_kN', 'Speed_km/h', 'Spacing_m', ...
        'Modulus_of_Elasticity_GPa', 'Moment_of_Inertia_m4', 'Damping_Ratio'});
    % Append to the full dataset
    random_data = [random_data; vehicle_data];
end
% Define the filename for the Excel file
filename = 'vehicle_data7.xlsx';
% Write the random data to an Excel file
writetable(random_data, filename);
% Display a message
disp(['Excel file "', filename, '" with random data has been created successfully.']);

[Walter Roberson]

Your code defines probability of failure as failure_count / num_simulations . failure_count is not known until all of the entries in the data are analyzed. failure_count is inherently a scalar, so probability_of_failure is inherently a scalar.

There are three ways to plot a scalar against something else:

1. plot() the scalar against another scalar, specifying a marker type. The end result will be only a single dot on the graph
2. and 3. Use xline() or yline() to plot the constant scalar probability against the whole width or height of the graph. This graph is mostly useless.

Plotting the scalar probability of failure is mostly uninteresting.

[Walter Roberson]

% Load the data
filename = 'vehicle_data7.xlsx';  % Data file
data = readtable(filename);       % Load the data into a table

% Beam properties
L = 40;                    % Beam length (m)
h = 2;                     % Beam height (m)
allowable_stress = 240e6;   % Allowable stress in Pa (240 MPa)
g = 9.81;                  % Gravitational acceleration (m/s^2)

% Constants
rho_steel = 7850;           % Density of steel (kg/m^3)

% Number of simulations (equal to number of rows in dataset)
num_simulations = height(data);

% Initialize counters for failures
failure_count = 0;
stresses = zeros(1, num_of_simulations);

% Monte Carlo simulation loop
for sim = 1:num_simulations
    % Extract vehicle data
    vehicle_type = data.VehicleType{sim}; % Vehicle type (e.g., sedan, bus, truck)
    weight = data.Weight_kN(sim) * 1e3;   % Convert from kN to N
    speed = data.Speed_km_h(sim) / 3.6;   % Convert from km/h to m/s
    spacing = data.Spacing_m(sim);        % Spacing (m)
    E_veh = data.Modulus_of_Elasticity_GPa(sim) * 1e9; % Convert from GPa to Pa
    I_veh = data.Moment_of_Inertia_m4(sim); % Moment of inertia (m^4)
    damping_ratio = data.Damping_Ratio(sim); % Damping ratio from dataset

    % Beam's cross-sectional area and mass per unit length
    b = I_veh * 12 / h^3;                 % Calculate beam width from moment of inertia
    A = b * h;                            % Cross-sectional area (m^2)
    m = rho_steel * A;                    % Mass per unit length (kg/m)

    % Modal mass (corrected formula)
    m_n = (m * L) / 2;                    % Modal mass (kg)

    % Natural frequency for the first mode (omega_1)
    omega_1 = sqrt((E_veh * I_veh) / (m_n * L^4));  % Natural frequency for first mode (rad/s)

    % Stiffness for the first mode
    k_1 = m_n * omega_1^2;               % Modal stiffness (N/m)

    % Critical velocity (m/s)
    v_c = omega_1 * (L / pi);            % Critical velocity

    % Dynamic Amplification Factor (DAF) calculation
    DAF = 1 / sqrt((1 - (speed / v_c)^2)^2 + (2 * damping_ratio * (speed / v_c))^2);


    % Mode shape function for the first mode (phi_1)
    phi_1 = @(x) sin(pi * x / L);        % Mode shape for first mode (phi_1)
    F_1 = weight * phi_1(L / 2);         % Modal force due to moving load (N)

    % Maximum displacement
    w_max = (F_1 / k_1) * DAF;          % Maximum dynamic displacement (m)

    % Maximum stress (sigma_max)
    phi_1_dd = @(x) -(pi^2 / L^2) * sin(pi * x / L); % Second derivative of phi_1
    z = h / 2;                           % Distance from neutral axis to outer fiber (m)
    M_max = -E_veh * I_veh * phi_1_dd(L / 2) * w_max; % Maximum bending moment (N.m)
    sigma_max = M_max * z / I_veh;       % Maximum dynamic stress (Pa)
    stresses(sim) = sigma_max;
    % Compare sigma_max with allowable stress to determine failure probability
    if sigma_max > allowable_stress
        failure_count = failure_count + 1;  % Increment failure counter
    end
end

% Calculate failure probability
failure_probability = failure_count / num_simulations;

% Display results
fprintf('Number of Simulations: %d\n', num_simulations);
fprintf('Number of Failures: %d\n', failure_count);
fprintf('Failure Probability: %.4f\n', failure_probability);

plot(stresses, failure_probability, '*-')

The above code plots failure probability versus stresses. You will notice that the plot is mostly useless because failure probability is a scalar.