% RBDO: Berm Breakwater with Multi-Objective Optimization
% 1. INITIALIZE UQLAB
clearvars
rng(1, 'twister')
uqlab

% 2. COMPUTATIONAL MODEL
% 2.3 Limit state function (constraint)
ModelOpts_Constraint.mFile = 'uq_Ex08_BermBreakwater_constraint';
myModel = uq_createModel(ModelOpts_Constraint);
RBDOOpts.LimitState.Model = myModel; % Correct field: Model, not mFile

% 2.1 Define the cost model (single output: Cost)
ModelOpts_Cost.mFile = 'uq_Ex08_BermBreakwater_cost';
myCostModel = uq_createModel(ModelOpts_Cost);
RBDOOpts.Cost.Model = myCostModel;

% 2.2 Define the objectives model for multi-objective optimization
ModelOpts_Objectives.mFile = 'uq_Ex08_BermBreakwater_objectives';
ModelOpts_Objectives.OutputDim = 3; % Three objectives: Cost, Carbon, Robustness
ModelOpts_Objectives.Parameters.ConstraintModel = myModel; % Now 'myModel' is defined
ModelOpts_Objectives.Parameters.EnvInput = myInput; % Pass environmental input
myObjectivesModel = uq_createModel(ModelOpts_Objectives);
RBDOOpts.Objectives = myObjectivesModel;


% 3. PROBABILISTIC INPUT MODEL
% 3.1 Design variables
RBDOOpts.Optim.Bounds = [14.3 0.9 2.6 1.1; 23.1 5 2.9 6];
RBDOOpts.Optim.StartingPoint = [17 3 2.75 2.83];

B_mean = RBDOOpts.Optim.StartingPoint(1);
hb_mean = RBDOOpts.Optim.StartingPoint(2);
Ns_mean = RBDOOpts.Optim.StartingPoint(3);
hf_mean = RBDOOpts.Optim.StartingPoint(4);

RBDOOpts.Input.DesVar(1).Name = '$B[m]$';
RBDOOpts.Input.DesVar(1).Type = 'Gaussian';
RBDOOpts.Input.DesVar(1).Moments = [B_mean, B_mean * 0.05];

RBDOOpts.Input.DesVar(2).Name = '$h_{b}[m]$';
RBDOOpts.Input.DesVar(2).Type = 'Gaussian';
RBDOOpts.Input.DesVar(2).Moments = [hb_mean, hb_mean * 0.05];

RBDOOpts.Input.DesVar(3).Name = '$N_{s}[-]$';
RBDOOpts.Input.DesVar(3).Type = 'Gaussian';
RBDOOpts.Input.DesVar(3).Moments = [Ns_mean, Ns_mean * 0.05];

RBDOOpts.Input.DesVar(4).Name = '$h_{f}[m]$';
RBDOOpts.Input.DesVar(4).Type = 'Gaussian';
RBDOOpts.Input.DesVar(4).Moments = [hf_mean, hf_mean * 0.05];

% 3.2 Environmental variables
InputOpts.Marginals(1).Name = 'Hs';
InputOpts.Marginals(1).Type = 'Lognormal';
InputOpts.Marginals(1).Moments = [5.03 0.7545];

InputOpts.Marginals(2).Name = 'Tp';
InputOpts.Marginals(2).Type = 'Lognormal';
InputOpts.Marginals(2).Moments = [13.6 2.04];

InputOpts.Marginals(3).Name = 'ro_s';
InputOpts.Marginals(3).Type = 'Gaussian';
InputOpts.Marginals(3).Moments = [2.45 0.049];

InputOpts.Marginals(4).Name = 'ro_w';
InputOpts.Marginals(4).Type = 'Gaussian';
InputOpts.Marginals(4).Moments = [1.025 0.01025];

InputOpts.Marginals(5).Name = 'd';
InputOpts.Marginals(5).Type = 'Gaussian';
InputOpts.Marginals(5).Moments = [16.3 0.518];

InputOpts.Marginals(6).Name = 'N';
InputOpts.Marginals(6).Type = 'Gaussian';
InputOpts.Marginals(6).Moments = [2978 992];

InputOpts.Marginals(7).Name = 'fg';
InputOpts.Marginals(7).Type = 'Gaussian';
InputOpts.Marginals(7).Moments = [1.5 0.075];

InputOpts.Marginals(8).Name = 'DD';
InputOpts.Marginals(8).Type = 'Gaussian';
InputOpts.Marginals(8).Moments = [1.0 0.05];

InputOpts.Marginals(9).Name = 'CR';
InputOpts.Marginals(9).Type = 'Gaussian';
InputOpts.Marginals(9).Moments = [1.0 0.1];

InputOpts.Marginals(10).Name = 'CO1';
InputOpts.Marginals(10).Type = 'Gaussian';
InputOpts.Marginals(10).Moments = [0.09 0.0135];

InputOpts.Marginals(11).Name = 'CO2';
InputOpts.Marginals(11).Type = 'Gaussian';
InputOpts.Marginals(11).Moments = [1.5 0.15];

% 3.3 Copula (disabled for testing to avoid Inference issues)
% InputOpts.Copula(1) = uq_VineCopula(...
%     'DVine', [1 2 3], ...
%     {'Gumbel', 'Clayton', 'Independent'}, {1.2467, 0.3613, []}, [0 180 0]);
% InputOpts.Copula(1).Variables = [1 2 5];

% Create input
myInput = uq_createInput(InputOpts);
RBDOOpts.Input.EnvVar = myInput;

% 4. OPTIMIZATION SETUP
% 4.1 Target Failure Probability (realistic value)
RBDOOpts.TargetPf = [0.01 0.01];

% 5. MULTI-OBJECTIVE RBDO SETUP
RBDOOpts.Type = 'RBDO';
RBDOOpts.Method = 'two-level';
RBDOOpts.MultiObjective = true;
RBDOOpts.Optim.Method = 'gamultiobj';
RBDOOpts.Optim.gamultiobj.PopulationSize = 200;
RBDOOpts.Optim.gamultiobj.MaxGenerations = 200;
RBDOOpts.Optim.gamultiobj.ParetoFraction = 0.7;
% RBDOOpts.Optim.gamultiobj.FunctionTolerance = 1e-6;
% RBDOOpts.Optim.gamultiobj.UseParallel = true;  % ??? ????? ??? ?? Parallel ???????? ???

% 6. RELIABILITY ANALYSIS SETUP
RBDOOpts.Reliability.Simulation.MaxSampleSize = 5e4;

% 7. RUN THE MULTI-OBJECTIVE RBDO ANALYSIS
myRBDO_MO = uq_createAnalysis(RBDOOpts);

% 8. PRINT AND DISPLAY RESULTS
uq_print(myRBDO_MO)
uq_display(myRBDO_MO)

% 9. POST-PROCESS PARETO FRONT
ParetoSolutions = myRBDO_MO.Results.ParetoSet;
ParetoObjectives = myRBDO_MO.Results.ParetoObjectives;

fprintf('Number of Pareto-optimal solutions: %d\n', size(ParetoSolutions, 1));
for i = 1:size(ParetoSolutions, 1)
    fprintf('Solution %d:\n', i);
    fprintf('  Design Variables [B, hb, Ns, hf]: [%.2f, %.2f, %.2f, %.2f]\n', ParetoSolutions(i,:));
    fprintf('  Objectives [Cost, Carbon, Robustness]: [%.2f, %.2f, %.2f]\n', ParetoObjectives(i,:));
end

% Plot Pareto front (3D scatter and 2D projections)
figure;
subplot(2, 2, 1);
scatter3(ParetoObjectives(:,1), ParetoObjectives(:,2), ParetoObjectives(:,3), 'filled');
xlabel('Lifecycle Cost'); ylabel('Carbon Footprint (kg CO_2)'); zlabel('Robustness (1/CoV)');
title('Pareto Front (3D)');
grid on;

subplot(2, 2, 2);
scatter(ParetoObjectives(:,1), ParetoObjectives(:,2), 'filled');
xlabel('Lifecycle Cost'); ylabel('Carbon Footprint');
title('Cost vs Carbon');
grid on;

subplot(2, 2, 3);
scatter(ParetoObjectives(:,1), ParetoObjectives(:,3), 'filled');
xlabel('Lifecycle Cost'); ylabel('Robustness');
title('Cost vs Robustness');
grid on;

subplot(2, 2, 4);
scatter(ParetoObjectives(:,2), ParetoObjectives(:,3), 'filled');
xlabel('Carbon Footprint'); ylabel('Robustness');
title('Carbon vs Robustness');
grid on;

% 10. CALCULATE AND DISPLAY FINAL DESIGN METRICS
Xstar = ParetoSolutions(1,:);
Nso = Xstar(3);
hbo = Xstar(2);
hfo = Xstar(4);
Dno = 5.03 / (((2.45/1.025) - 1) * Nso);
HSTp = (Nso * sqrt(13.6 * sqrt(9.81/Dno)));
Nd = 9 * 3600 / 0.8 / 13.6;

REC = (((0.02067 * (HSTp^1.971) * ((16.3/Dno)^-0.1087) * log((Nd*16.3)/(3000*Dno))) + ...
        (21133 * (HSTp^2.094) * (hbo/5.03) * (16.3/Dno) * (log(1.5)^(7.187^1.5))) - ...
        (81.74 * (10^-5) * (1.378^(hbo/5.03)) * (HSTp^3) * (Nd/3000) * ((hbo/5.03)^2)) - 1.245) * Dno);
Rc = hfo + 17.8 + 2 * Dno - 16.3;
M = 2.45 * (Dno^3);

fprintf('Selected Pareto Solution:\n');
fprintf('The value of Dno (m) is: %.2f\n', Dno);
fprintf('The value of M (ton) is: %.2f\n', M);
fprintf('The value of Rec (m) is: %.2f\n', REC);
fprintf('The value of Rc (m) is: %.2f\n', Rc);

% 11. CHECK CONSTRAINTS
nSamples = 5e4;
envSamples = uq_getSample(myInput, nSamples);
X_full = [repmat(Xstar, nSamples, 1), envSamples];
limitStateOutput = uq_evalModel(myModel, X_full);
Pf_REC = mean(limitStateOutput(:,1) <= 0);
Pf_Overtopping = mean(limitStateOutput(:,2) <= 0);
fprintf('Failure Probabilities for Selected Solution:\n');
fprintf('  REC: %.4f (Target: %.4f)\n', Pf_REC, RBDOOpts.TargetPf(1));
fprintf('  Overtopping: %.4f (Target: %.4f)\n', Pf_Overtopping, RBDOOpts.TargetPf(2));