Skip to content

Commit

Permalink
Merge pull request #47 from symbiotic-engineering/iris-design-of-expe…
Browse files Browse the repository at this point in the history
…riments

Iris design of experiments
  • Loading branch information
rebeccamccabe authored May 15, 2024
2 parents c10a49c + 2783276 commit 6cc51a0
Show file tree
Hide file tree
Showing 5 changed files with 160 additions and 143 deletions.
5 changes: 4 additions & 1 deletion .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -43,4 +43,7 @@ code-coverage/
test-results/

# mat data files
*.mat
*.mat

# excel cache
~$*.xlsx
131 changes: 0 additions & 131 deletions dev/experiments.m

This file was deleted.

13 changes: 11 additions & 2 deletions mdocean/inputs/var_bounds.m
Original file line number Diff line number Diff line change
@@ -1,40 +1,48 @@
function b = var_bounds()

% outer diameter of float (m)
b.D_f_min = 6; % because D_s_nom = 6, a consequence of the spar natural frequency constraint
b.D_f_max = 40;
b.D_f_nom = 20;
b.D_f_start = 20;

% inner diameter to outer diameter of float ratio (-)
b.D_s_ratio_min = 0.01;
b.D_s_ratio_max = 0.99;
b.D_s_ratio_nom = 6/20;
b.D_s_ratio_start = 6/20;

% height of float to outer diameter of float ratio (-)
b.h_f_ratio_min = .1;
b.h_f_ratio_max = 10;
b.h_f_ratio_nom = 4/20;
b.h_f_ratio_start = 4/20;

% draft of spar to height of spar ratio (-)
b.T_s_ratio_min = 0.01;
b.T_s_ratio_max = 0.99;
b.T_s_ratio_nom = 35/44;
b.T_s_ratio_start = 35/44;

% maximum powertrain force (MN)
b.F_max_min = 0.01;
b.F_max_max = 100;
b.F_max_nom = 5;
b.F_max_start = 5;

% powertrain damping (MN / (m/s))
b.B_p_min = .1;
b.B_p_max = 50;
b.B_p_nom = 10;
b.B_p_start = 0.5;

% natural frequency (rad/s)
b.w_n_min = .01;%2*pi/p.T(find(any(p.JPD > 0),1,'last')); % min wave frequency that has any energy
b.w_n_max = 40;%2*pi/p.T(find(any(p.JPD > 0),1,'first')); % max wave frequency that has any energy
b.w_n_nom = 0.8;
b.w_n_start = 0.8;

% material index (-)
b.M_min = 1;
b.M_max = 3;
b.M_nom = 1;
Expand All @@ -50,9 +58,10 @@
'F_max',b.F_max_start,'D_int',b.B_p_start,'w_n',b.w_n_start);

b.var_names = {'D_f','D_s_ratio','h_f_ratio','T_s_ratio','F_max','B_p','w_n','M'};
b.var_names_pretty = {'D_f','D_s/D_f','h_f/D_f','T_s/h_s','F_{max}','B_p','\omega_n','M'};
b.constraint_names = {'float_too_heavy','float_too_light','spar_too_heavy','spar_too_light',...
'stability','FOS_float_yield','FOS_col_yield','FOS_plate','FOS_col_buckling',...
'pos_power','spar_damping','spar_height','LCOE_max','irrelevant_max_force'};
'stability','FOS_float_yield','FOS_col_yield','FOS_plate','FOS_col_buckling',...
'pos_power','spar_damping','spar_height','LCOE_max','irrelevant_max_force'};

% modify nominal control inputs to so power and force matches actual
[F_max_nom, B_p_nom, w_n_nom] = find_nominal_inputs(b, false);
Expand Down
18 changes: 9 additions & 9 deletions mdocean/optimization/multiobjective/pareto_curve_heuristics.m
Original file line number Diff line number Diff line change
Expand Up @@ -86,21 +86,21 @@ function pareto_curve_heuristics()
x_best_LCOE,x_best_Pvar,x_nom,x_balanced,idxo,showSingleObj,showImages,p)
figure
% overall pareto front
plot(LCOE(idxo),Pvar(idxo),'bs','MarkerFaceColor','b')
plot(LCOE(idxo),Pvar(idxo),'bs','MarkerFaceColor','b','HandleVisibility','off')
hold on

% utopia point
plot(minLCOE,minPvar,'gp','MarkerFaceColor','g','MarkerSize',20)
plot(minLCOE,minPvar,'gp','MarkerFaceColor','g','MarkerSize',20,'HandleVisibility','off')

% RM3 nominal reference
plot(LCOE_nom,P_var_nom,'rd')
plot(LCOE_nom_sim,P_var_nom_sim,'rs')
plot(LCOE_nom,P_var_nom,'rd','HandleVisibility','off')
plot(LCOE_nom_sim,P_var_nom_sim,'rs','HandleVisibility','off')

if showSingleObj
% black squares for 3 ref points
plot(minLCOE,Pvar(idx_best_LCOE),'ks')
plot(LCOE(idx_best_Pvar),minPvar,'ks')
plot(LCOE_balanced,P_var_balanced,'ks')
plot(minLCOE,Pvar(idx_best_LCOE),'ks','HandleVisibility','off')
plot(LCOE(idx_best_Pvar),minPvar,'ks','HandleVisibility','off')
plot(LCOE_balanced,P_var_balanced,'ks','HandleVisibility','off')
end

% axis labels
Expand All @@ -112,7 +112,7 @@ function pareto_curve_heuristics()
improvePlot

% solar reference
plot(LCOE_solar, P_var_solar,'o','MarkerSize',12,'MarkerEdgeColor',[1, .87, .2],'MarkerFaceColor',[1, .87, .2])
plot(LCOE_solar, P_var_solar,'o','MarkerSize',12,'MarkerEdgeColor',[1, .87, .2],'MarkerFaceColor',[1, .87, .2],'HandleVisibility','off')
% for the yellow color to work, do not use improvePlot below here

% text labels
Expand All @@ -128,7 +128,7 @@ function pareto_curve_heuristics()
text(LCOE(idx_best_Pvar)+.03,Pvar(idx_best_Pvar)-3,'Least Variable','FontSize',sz)
text(LCOE_balanced-.15,P_var_balanced+5,'Balanced Design','FontSize',sz)
end

if showImages
mini_plot_size = [.2 .22];
% small corner pictures of best geometries
Expand Down
136 changes: 136 additions & 0 deletions mdocean/simulation/run/experiments.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@

% Runs one-at-a-time design of experiments

clear;clc;close all

p = parameters();
b = var_bounds();

DOE_strategy = 'ratios'; % 'sample' or 'bounds' or 'ratios'
n = 20;

if strcmp(DOE_strategy,'sample')
X = [ 20 10 30; % outer diameter of float
.3 .1 .5; % inner diameter ratio of float
1.5 1.2 2; % outer diameter of reaction plate
1 2 3; % material
10 1 20; % Number of WECs in array
10 5 50 % D_int
2*pi/7 2*pi/8 2*pi/9
];
elseif strcmp(DOE_strategy,'bounds')
X = [ b.D_f_nom, linspace(b.D_f_min, b.D_f_max, n);
b.D_s_ratio_nom, linspace(b.D_s_ratio_min, b.D_s_ratio_max, n);
b.h_f_ratio_nom, linspace(b.h_f_ratio_min, b.h_f_ratio_max, n);
b.T_s_ratio_nom, linspace(b.T_s_ratio_min, b.T_s_ratio_max, n);
b.F_max_nom, linspace(b.F_max_min, b.F_max_max, n);
b.B_p_nom, linspace(b.B_p_min, b.B_p_max, n);
b.w_n_nom, linspace(b.w_n_min, b.w_n_max, n) ];
ratios = X./X(:,1);
elseif strcmp(DOE_strategy,'ratios')
ratios = logspace(log10(1/3),log10(3),n);
ratios = [1, ratios(ratios~=1)];
X = [b.D_f_nom b.D_s_ratio_nom b.h_f_ratio_nom b.T_s_ratio_nom b.F_max_nom b.B_p_nom b.w_n_nom b.M_nom]' * ratios;
end

X_nom = X(:,1);
design_size = size(X);
num_vars = design_size(1);
num_vals_per_var = design_size(2);
num_vals_swept = num_vals_per_var - 1; % don't sweep 1st value of each variable (corresponding to ratio 1)

% don't sweep material
idx_material = 8;
X(idx_material,:) = b.M_nom * ones(1,n+1);
num_vars_swept = num_vars - 1;

% initialize variables
LCOE = X*inf;
P_var = X*inf;
opt_idx = zeros(num_vars,1);
recommended = zeros(num_vars,2);
number_runs = 1 + num_vars_swept * num_vals_swept; % nominal run plus all sweeps
failed = cell(number_runs,1);
power = zeros(number_runs,1);
FOS = zeros(number_runs,1);
X_ins = zeros(number_runs, num_vars);
design = 0;

% run design of experiments
for i = 1:num_vars_swept
X_in = X_nom;
for j = 1:num_vals_per_var
if i == 1 || j~=1 % prevent rerunning nominal multiple times
changed_entry = X(i,j);
if ~isnan(changed_entry)
design = design+1;
X_in(i) = changed_entry;
X_ins(design,:) = X_in;
[LCOE_temp, P_var_temp, P_matrix, g, val] = simulation(X_in,p);

% only add to results if first 12 constraints are feasible
g(13) = 0;
g(14) = 0;
[feasible, which_failed] = is_feasible(g, b);
if feasible
LCOE(i,j) = LCOE_temp;
P_var(i,j) = P_var_temp;
else
LCOE(i,j) = NaN;
P_var(i,j) = NaN;
end
failed{design} = which_failed;
end
end
end
[~, opt_idx(i)] = min(LCOE(i,:));
recommended(i,:) = [X(i,opt_idx(i)), opt_idx(i)];
end

% create table for display
results = array2table(X_ins, 'VariableNames', b.var_names);
LCOE = LCOE';
P_var = P_var';
results = addvars(results, round(LCOE(LCOE~=Inf),2), round(P_var(P_var~=Inf),1), failed, ...
'NewVariableNames', {'LCOE ($/kWh)','c_v (%)','Failed Constraints'});
disp(results)

% plot pareto curve for comparison
pareto_curve_heuristics()
figure(2)
plot(LCOE, P_var, '*--')
xlabel('LCOE')
ylabel('P_{var}')
title('Design of Experiments Pareto Front')
l = legend(b.var_names_pretty);
improvePlot
l.Location = 'bestoutside';

%% sensitivities plot
[ratios_sorted,idx] = sort(ratios);
LCOE(1,:) = LCOE(1,1); % fill in nominal LCOE results for each DV where it wasn't repeatedly tested
P_var(1,:) = P_var(1,1);

figure
t = tiledlayout(2,1);
t.TileSpacing = 'compact';

ax1 = nexttile(1);
plot(ratios_sorted,LCOE(idx,:).')
ylabel('LCOE ($/kWh)')
grid on
ax2 = nexttile(2);
plot(ratios_sorted,P_var(idx,:).')
ylabel('Power c_v (%)')
grid on

title(t,'Design of Experiments Results','FontWeight','bold','FontSize',20)
l = legend(b.var_names_pretty);
l.Location = 'bestoutside';
grid on
linkaxes([ax1,ax2],'x');
xlabel('Design Variable Ratio (-)')
xticklabels(ax1,{})
xticks(ax2,xticks(ax1))
improvePlot

0 comments on commit 6cc51a0

Please sign in to comment.