%% ROVER WHEEL CONFIGURATION ANALYSIS
% Reads the FEM simulation table (Patran/Nastran),
% builds the worst-case summary for each design,
% filters feasible configurations,
% computes weighted ranking,
% identifies the Pareto front,
% selects 3 significant cases.
%
% NOTES:
% 1) Save your table in Excel.
% 2) Modify fileName below.
% 3) Check the design limits (sigma_allow, u_lim, etc.).


clear; clc; close all;

%% ===================== USER INPUT =====================
fileName = 'TAB_tweel_final.xlsx';   
sheetName = 1;

% Analysis mode:
% "design_worst_case" = recommended if the same geometry appears under multiple loads
% "row"               = uses each row as an independent configuration
analysisMode = "row";

% ---- Design constraints ----
sigma_allow = 264e6;       % [Pa] allowable material stress
SF_required = 1.50;        % [-] required safety factor
u_lim = 0.030;             % [m] maximum allowable vertical deflection
mass_max = inf;            % [kg] set a number if you want a mass constraint
stiffness_min = 0;         % units consistent with your column

% Durability index:
% true  -> higher = better, enforce durability >= durability_limit
% false -> lower = better, enforce durability <= durability_limit
paramsDurability.higherIsBetter = true;
paramsDurability.limit = 0.0;

% Weights for the final ranking (must sum to 1)
weights.stress     = 0.35;
weights.deflection = 0.25;
weights.mass       = 0.20;
weights.stiffness  = 0.10;
weights.durability = 0.10;

outPrefix = 'risultati_ruote';

%% 1) Read table
opts = detectImportOptions(fileName, 'Sheet', sheetName, 'VariableNamingRule', 'preserve');
Traw = readtable(fileName, opts);

fprintf('\nTabella letta: %d righe, %d colonne\n', height(Traw), width(Traw));
disp('Colonne trovate nel file:');
disp(string(Traw.Properties.VariableNames)');

%% 2) Mapping of your table columns
vars.rExt            = pickVarName(Traw, {'R ext', 'Rext', 'R esterno'});
vars.rInt            = pickVarName(Traw, {'R int', 'Rint', 'R interno'});
vars.width           = pickVarName(Traw, {'Width', 'Larghezza'});
vars.tExt            = pickVarName(Traw, {'t ext', 't_ext', 'spessore esterno'});
vars.tInt            = pickVarName(Traw, {'t int', 't_int', 'spessore interno'});
vars.tRinf           = pickVarName(Traw, {'t rinf', 't_rinf', 'spessore rinforzo'});
vars.vehicleWeight   = pickVarName(Traw, {'Vehicle weight', 'Peso veicolo'});
vars.load            = pickVarName(Traw, {'Load', 'Carico'});
vars.nRinf           = pickVarName(Traw, {'N° rinf', 'N rinf', 'Nrinf', 'numero rinforzi'});
vars.matExtInt       = pickVarName(Traw, {'Material ext & int', 'Material ext int', 'materiale est int'});
vars.matRinf         = pickVarName(Traw, {'Material rinf', 'materiale rinforzo'});
vars.verticalDefl    = pickVarName(Traw, {'Vertical deflection', 'Deflessione verticale'});
vars.maximumStress   = pickVarName(Traw, {'Maximum stress', 'Stress max', 'Tensione massima'});
vars.wheelMass       = pickVarName(Traw, {'Wheel mass', 'Massa ruota'});
vars.wheelStiffness  = pickVarName(Traw, {'Wheel stiffness', 'Rigidezza ruota'});
vars.durabilityIndex = pickVarName(Traw, {'Durability index', 'Indice di durabilita', 'Durability'});

%% 3) Build cleaned table
T = table();
T.Rext               = getColumn(Traw, vars.rExt);
T.Rint               = getColumn(Traw, vars.rInt);
T.Width              = getColumn(Traw, vars.width);
T.tExt               = getColumn(Traw, vars.tExt);
T.tInt               = getColumn(Traw, vars.tInt);
T.tRinf              = getColumn(Traw, vars.tRinf);
T.VehicleWeight      = getColumn(Traw, vars.vehicleWeight);
T.Load               = getColumn(Traw, vars.load);
T.NRinf              = getColumn(Traw, vars.nRinf);
T.MatExtInt          = getColumn(Traw, vars.matExtInt);
T.MatRinf            = getColumn(Traw, vars.matRinf);
T.VerticalDeflection = getColumn(Traw, vars.verticalDefl);
T.MaximumStress      = getColumn(Traw, vars.maximumStress);
T.WheelMass          = getColumn(Traw, vars.wheelMass);
T.WheelStiffness     = getColumn(Traw, vars.wheelStiffness);
T.DurabilityIndex    = getColumn(Traw, vars.durabilityIndex);

keyMask = ~isnan(asNumeric(T.MaximumStress)) & ~isnan(asNumeric(T.VerticalDeflection));
T = T(keyMask,:);

tol = 1e-12;  % numerical tolerance
maskWrongThickness = ...
    abs(T.tInt  - 0.003) < tol | ...
    abs(T.tExt  - 0.003) < tol;

fprintf('Configurazioni eliminate per spessore = 0.003 m: %d\n', sum(maskWrongThickness));

T = T(~maskWrongThickness,:);

fprintf('\nDopo pulizia: %d righe\n', height(T));

%% 4) Row-wise analysis or worst-case by design
if analysisMode == "row"
    A = T;
    A.DesignID = (1:height(A))';
    A.nLoadCases = ones(height(A),1);

    stressCol = 'MaximumStress';
    deflCol   = 'VerticalDeflection';
    massCol   = 'WheelMass';
    stiffCol  = 'WheelStiffness';
    durCol    = 'DurabilityIndex';

elseif analysisMode == "design_worst_case"
    designVars = {'Rext','Rint','Width','tExt','tInt','tRinf','NRinf','MatExtInt','MatRinf'};
    [designUnique, ~, G] = unique(T(:, designVars), 'rows', 'stable');

    A = designUnique;
    A.DesignID = (1:height(designUnique))';
    A.nLoadCases      = splitapply(@numel, T.MaximumStress, G);
    A.MaxVehicleWeight = splitapply(@maxnum, asNumeric(T.VehicleWeight), G);
    A.MaxLoad          = splitapply(@maxnum, asNumeric(T.Load), G);
    A.WorstStress      = splitapply(@maxnum, asNumeric(T.MaximumStress), G);
    A.WorstDeflection  = splitapply(@maxnum, asNumeric(T.VerticalDeflection), G);
    A.Mass             = splitapply(@meannum, asNumeric(T.WheelMass), G);
    A.MinStiffness     = splitapply(@minnum, asNumeric(T.WheelStiffness), G);

    if paramsDurability.higherIsBetter
        A.WorstDurability = splitapply(@minnum, asNumeric(T.DurabilityIndex), G);
    else
        A.WorstDurability = splitapply(@maxnum, asNumeric(T.DurabilityIndex), G);
    end

    stressCol = 'WorstStress';
    deflCol   = 'WorstDeflection';
    massCol   = 'Mass';
    stiffCol  = 'MinStiffness';
    durCol    = 'WorstDurability';
else
    error('analysisMode non valido. Usa "row" oppure "design_worst_case".');
end

%% 5) Constraints
A.SigmaLimit = repmat(sigma_allow / SF_required, height(A), 1);
A.StressRatio = asNumeric(A.(stressCol)) ./ A.SigmaLimit;
A.DeflectionRatio = asNumeric(A.(deflCol)) ./ u_lim;
A.StiffnessRatio = stiffness_min ./ max(asNumeric(A.(stiffCol)), eps);
A.MassRatio = asNumeric(A.(massCol)) ./ mass_max;

A.StressOK = asNumeric(A.(stressCol)) <= A.SigmaLimit;
A.DeflectionOK = asNumeric(A.(deflCol)) <= u_lim;
A.MassOK = asNumeric(A.(massCol)) <= mass_max;
A.StiffnessOK = asNumeric(A.(stiffCol)) >= stiffness_min;

if paramsDurability.higherIsBetter
    A.DurabilityOK = asNumeric(A.(durCol)) >= paramsDurability.limit;
    A.DurabilityMargin = asNumeric(A.(durCol)) - paramsDurability.limit;
else
    A.DurabilityOK = asNumeric(A.(durCol)) <= paramsDurability.limit;
    A.DurabilityMargin = paramsDurability.limit - asNumeric(A.(durCol));
end

A.Feasible = A.StressOK & A.DeflectionOK & A.MassOK & A.StiffnessOK & A.DurabilityOK;

fprintf('\nConfigurazioni/design totali: %d\n', height(A));
fprintf('Configurazioni/design ammissibili: %d\n', sum(A.Feasible));

if ~any(A.Feasible)
    error(['Nessuna configurazione e ammissibile con i vincoli attuali. ' ...
           'Rivedi limiti/pesi oppure controlla la tabella.']);
end

%% 6) Weighted ranking on feasible configurations only
Feas = A(A.Feasible,:);

nStress = normalizeCost(asNumeric(Feas.(stressCol)));
nDefl   = normalizeCost(asNumeric(Feas.(deflCol)));
nMass   = normalizeCost(asNumeric(Feas.(massCol)));
nStiff  = normalizeBenefit(asNumeric(Feas.(stiffCol)));
if paramsDurability.higherIsBetter
    nDur = normalizeBenefit(asNumeric(Feas.(durCol)));
else
    nDur = normalizeCost(asNumeric(Feas.(durCol)));
end

Feas.Score = weights.stress*nStress + ...
             weights.deflection*nDefl + ...
             weights.mass*nMass + ...
             weights.stiffness*nStiff + ...
             weights.durability*nDur;

%% 7) Pareto front
obj = [asNumeric(Feas.(stressCol)), ...
       asNumeric(Feas.(deflCol)), ...
       asNumeric(Feas.(massCol)), ...
       -asNumeric(Feas.(stiffCol))];
if paramsDurability.higherIsBetter
    obj = [obj, -asNumeric(Feas.(durCol))];
else
    obj = [obj, asNumeric(Feas.(durCol))];
end
Feas.Pareto = paretoNondominated(obj);

%% 8) Final ranking
primarySort = sortrows(Feas, {'Score', stressCol, deflCol, massCol}, {'ascend','ascend','ascend','ascend'});

A.Score = nan(height(A),1);
A.Pareto = false(height(A),1);
[tf, loc] = ismember(A.DesignID, Feas.DesignID);
A.Score(tf) = Feas.Score(loc(tf));
A.Pareto(tf) = Feas.Pareto(loc(tf));

best = primarySort(1,:);

%% 9) Selection of 3 significant cases
case1 = best;

case1.Criticality = max(case1.StressRatio, case1.DeflectionRatio);
if isfinite(paramsDurability.limit)
    critDur = 1 ./ max(case1.DurabilityMargin, 1e-12);
    critDur = critDur ./ max(critDur);
    case1.Criticality = max(case1.Criticality, critDur);
end

usedIDs = case1.DesignID;

Feas.Criticality = max([Feas.StressRatio, Feas.DeflectionRatio], [], 2);
if isfinite(paramsDurability.limit)
    critDur = 1 ./ max(Feas.DurabilityMargin, 1e-12);
    critDur = critDur / max(critDur);
    Feas.Criticality = max([Feas.Criticality, critDur], [], 2);
end
[~, idxCrit] = sort(Feas.Criticality, 'descend');
case2 = pickDistinctRow(Feas, idxCrit, usedIDs);
usedIDs = [usedIDs; case2.DesignID]; %#ok<AGROW>

[~, idxMassSort] = sort(asNumeric(Feas.(massCol)), 'ascend');
case3 = pickDistinctRow(Feas, idxMassSort, usedIDs);

cases3 = [case1; case2; case3];

%% 10) Print results to screen
fprintf('\n===== MIGLIORE CONFIGURAZIONE =====\n');
disp(best);

fprintf('\n===== 3 CASI SIGNIFICATIVI =====\n');
varsToShow = {'DesignID','Rext','Rint','Width','tExt','tInt','tRinf','NRinf', ...
              'MatExtInt','MatRinf', stressCol, deflCol, massCol, stiffCol, durCol, ...
              'Score','Pareto','Feasible'};
varsToShow = varsToShow(ismember(varsToShow, cases3.Properties.VariableNames));
disp(cases3(:, varsToShow));

%% 11) Plots
figure('Name','Massa vs Stress');
scatter(asNumeric(A.(massCol)), asNumeric(A.(stressCol)), 40, A.Feasible, 'filled');
grid on; xlabel('Massa [kg]'); ylabel('Stress [Pa]');
title('Massa vs Stress');
colormap(parula(2));
cb = colorbar; cb.Ticks = [0.25 0.75]; cb.TickLabels = {'Non feasible','Feasible'};

figure('Name','Pareto: Deflessione vs Stress');
feasDefl = asNumeric(Feas.(deflCol));
feasStress = asNumeric(Feas.(stressCol));
scatter(feasDefl, feasStress, 120, Feas.Score, 'filled'); hold on;
plot(feasDefl(Feas.Pareto), feasStress(Feas.Pareto), 'ko', 'MarkerSize', 15, 'LineWidth', 2);
grid on; xlabel('Deflection [m]',Interpreter='latex',FontSize=20); ylabel('Stress [Pa]',Interpreter='latex',FontSize=20);
% title('Feasible configurations and Pareto front');
cb = colorbar; cb.Label.String = 'Score';
legend('Feasible','Pareto','Location','southeast',Interpreter='latex',FontSize=25);
ax = gca;
ax.FontSize = 25;      % larger axis numbers
ax.LineWidth = 1;      % slightly thicker axes
ax.TickLabelInterpreter = "latex";

figure('Name','Spazio obiettivi');
scatter3(asNumeric(Feas.(massCol)), asNumeric(Feas.(deflCol)), asNumeric(Feas.(stressCol)), ...
         120, Feas.Score, 'filled');
grid on; xlabel('Massa [kg]',Interpreter='latex',FontSize=20); ylabel('Deflection [m]',Interpreter='latex',FontSize=20); zlabel('Stress [Pa]',Interpreter='latex',FontSize=20);
% title('Objective space of feasible configurations');
cb = colorbar; cb.Label.String = 'Score';
ax = gca;
ax.FontSize = 25;      % larger axis numbers
ax.LineWidth = 1;      % slightly thicker axes
ax.TickLabelInterpreter = "latex";



%% ===================== LOCAL FUNCTIONS =====================
function name = pickVarName(T, candidates)
    namesOriginal = string(T.Properties.VariableNames);
    namesNorm = normalizeNames(namesOriginal);
    candNorm = normalizeNames(string(candidates));
    for i = 1:numel(candNorm)
        idx = find(contains(namesNorm, candNorm(i)), 1, 'first');
        if ~isempty(idx)
            name = namesOriginal(idx);
            return;
        end
    end
    error('Non trovo una colonna compatibile con: %s', strjoin(string(candidates), ', '));
end

function out = normalizeNames(s)
    out = lower(regexprep(string(s), '[^a-zA-Z0-9]', ''));
end

function x = getColumn(T, varName)
    x = T.(varName);
    if iscell(x)
        try
            xnum = cellfun(@str2double, x);
            if ~all(isnan(xnum))
                x = xnum;
            end
        catch
        end
    elseif isstring(x)
        xnum = str2double(x);
        if ~all(isnan(xnum))
            x = xnum;
        end
    elseif iscategorical(x)
        xs = string(x);
        xnum = str2double(xs);
        if ~all(isnan(xnum))
            x = xnum;
        else
            x = cellstr(xs);
        end
    end
end

function x = asNumeric(x)
    if isnumeric(x)
        return;
    end
    if iscell(x)
        x = cellfun(@str2double, x);
    elseif isstring(x)
        x = str2double(x);
    elseif iscategorical(x)
        x = str2double(string(x));
    else
        x = double(x);
    end
end

function y = normalizeCost(x)
    x = asNumeric(x);
    xmin = min(x, [], 'omitnan');
    xmax = max(x, [], 'omitnan');
    if abs(xmax - xmin) < eps
        y = zeros(size(x));
    else
        y = (x - xmin) ./ (xmax - xmin);
    end
end

function y = normalizeBenefit(x)
    x = asNumeric(x);
    xmin = min(x, [], 'omitnan');
    xmax = max(x, [], 'omitnan');
    if abs(xmax - xmin) < eps
        y = zeros(size(x));
    else
        y = (xmax - x) ./ (xmax - xmin);
    end
end

function mask = paretoNondominated(F)
    n = size(F,1);
    mask = true(n,1);
    for i = 1:n
        if ~mask(i), continue; end
        for j = 1:n
            if i == j, continue; end
            if all(F(j,:) <= F(i,:)) && any(F(j,:) < F(i,:))
                mask(i) = false;
                break;
            end
        end
    end
end

function row = pickDistinctRow(T, idxOrder, usedIDs)
    row = [];
    for k = 1:numel(idxOrder)
        cand = T(idxOrder(k),:);
        if ~ismember(cand.DesignID, usedIDs)
            row = cand;
            return;
        end
    end
    row = T(idxOrder(1),:);
end

function y = maxnum(x)
    y = max(x, [], 'omitnan');
end

function y = minnum(x)
    y = min(x, [], 'omitnan');
end

function y = meannum(x)
    y = mean(x, 'omitnan');
end