Skip to content

Commit

Permalink
Merge pull request #63 from symbiotic-engineering/iris-design-vars
Browse files Browse the repository at this point in the history 
add more desgin variables to vector X
  • Loading branch information
@rebeccamccabe
rebeccamccabe authored Dec 2, 2024
2 parents 6bc7f17 + 99782e5 commit ee986b1
Show file tree
Hide file tree
Showing 17 changed files with 316 additions and 127 deletions.
Binary file added dev/float_frame_stress.mlx
View file
Binary file not shown.
33 changes: 20 additions & 13 deletions mdocean/inputs/parameters.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,10 +45,13 @@
table("h","h",{100},"site",true,"water depth (m)");
table("JPD","JPD",{jpd(2:end,2:end)},"site",false,...
"joint probability distribution of wave (%)");
table("Hs","H_s",{jpd(2:end,1)},"site",true,"wave height (m)");
table("Hs_struct","H_{s,struct}",{11.9},"site",true,"100 year wave height (m)");
table("T","T",{jpd(1,2:end)},"site",true,"wave period (s)");
table("T_struct","T_{struct}",{17.1},"site",true,"100 year wave period (s)");
table("Hs","H_s",{jpd(2:end,1)},"site",true,"significant wave height (m)");
table("T","T",{jpd(1,2:end)},"site",true,"wave energy period (s)");
table("Hs_struct","H_{s,struct}",{[5 7 9 11.22 9 7 5]*1.9*sqrt(2)},...
"site",true,"100 year significant individual wave height (m)");
table("T_struct","T_{struct}",{[5.57 8.76 12.18 17.26 21.09 24.92 31.70]},...
"site",true,"100 year wave peak period (s)");
...% 100 year extreme heights and periods from Berg 2011
...
...% Materials: [ Structural Steel ASTM-A36 (Ductile)
...% Concrete (Brittle)
Expand All @@ -65,16 +68,16 @@
table("nu","\nu",{[0.36 0 0.29]},"structures",true,"Poisson's ratio (-)");
...
...% Thicknesses and structures
table("t_f_t","t_{f,t}",{0.50 * in2m},"structures",true,"float top thickness (m)");
table("t_f_r","t_{f,r}",{0.44 * in2m},"structures",true,"float radial wall thickness (m)");
table("t_f_c","t_{f,c}",{0.44 * in2m},"structures",true,...
"float circumferential gusset thickness (m)");
table("t_f_b","t_{f,b}",{0.56 * in2m},"structures",true,"float bottom thickness (m)");
table("t_s_r","t_{s,r}",{1.00 * in2m},"structures",true,"vertical column thickness (m)");
...% table("t_f_t","t_{f,t}",{0.50 * in2m},"structures",true,"float top thickness (m)");
...% table("t_f_r","t_{f,r}",{0.44 * in2m},"structures",true,"float radial wall thickness (m)");
...% table("t_f_c","t_{f,c}",{0.44 * in2m},"structures",true,...
...% "float circumferential gusset thickness (m)");
...% table("t_f_b","t_{f,b}",{0.56 * in2m},"structures",true,"float bottom thickness (m)");
...% table("t_s_r","t_{s,r}",{1.00 * in2m},"structures",true,"vertical column thickness (m)");
table("t_d_max","t_{d,max}",{1.00 * in2m},"structures",true,...
"max thickness of damping plate before making it hollow (m)");
table("t_d_tu","t_{d,tu}",{1.00 * in2m},"structures",true,...
"damping plate support tube radial wall thickness (m)");
...% table("t_d_tu","t_{d,tu}",{1.00 * in2m},"structures",true,...
...% "damping plate support tube radial wall thickness (m)");
table("D_d_tu","D_{d,tu}",{48.00 * in2m},"structures",true,...
"damping plate support tube diameter (m)");
table("theta_d_tu","\theta_{d,tu}",{atan(17.5/15)},"structures",true,...
Expand Down Expand Up @@ -143,7 +146,11 @@
table("spar_excitation_coeffs","spar excitation coeffs",{spar_exc},...
"dynamics",false,"spar excitation hydro coeffs from WAMIT for nominal RM3");
table("hydro","hydro",{readWAMIT(struct(),"rm3.out",[])},"dynamics",...
false,"function from WECSim")];
false,"function from WECSim");
table("F_heave_mult","F_{heave,mult}",{1.925*0.857},"dynamics",true,...
"multiplier to make heave force match with validation (-)")];
% 1.925 is required to make WAMIT match tank test, and 0.857 is
% required to make MEEM match WAMIT

T.Properties.VariableNames = cols;

Expand Down
84 changes: 73 additions & 11 deletions mdocean/inputs/var_bounds.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,10 @@
mode = '';
end

b.var_names = {'D_s','D_f','T_f_2','h_s','F_max','B_p','h_fs_clear','M'};
b.var_names_pretty = {'D_s','D_f','T_{f,2}','h_s','F_{max}','B_p','h_{fs,clear}','M'};
b.var_names = {'D_s','D_f','T_f_2','h_s','F_max','B_p','h_fs_clear',...
't_ft','t_fr','t_fc','t_fb','t_sr','t_dt','P_max','M'};
b.var_names_pretty = {'D_s','D_f','T_{f,2}','h_s','F_{max}','B_p','h_{fs,clear}',...
't_{ft}','t_{fr}','t_{fc}','t_{fb}','t_{sr}','t_{dt}','P_{max}','M'};

% inner diameter of float (m)
b.D_s_min = 0;
Expand All @@ -18,7 +20,7 @@

% outer diameter of float (m)
b.D_f_min = 1;
b.D_f_max = 40;
b.D_f_max = 50;
b.D_f_wecsim = 20;
b.D_f_nom = 20;
b.D_f_start = 20;
Expand Down Expand Up @@ -58,27 +60,85 @@
b.h_fs_clear_nom = 4; % p169 11/26/24
b.h_fs_clear_start = 4;

in2mm = 25.4;
% material thickness of float top (mm)
b.t_ft_min = 0.1 * in2mm;
b.t_ft_max = 1.0 * in2mm;
b.t_ft_nom = 0.5 * in2mm;
b.t_ft_wecsim = 0.5 * in2mm;
b.t_ft_start = 0.5 * in2mm;

% material thickness of float bottom (mm)
b.t_fr_min = 0.1 * in2mm;
b.t_fr_max = 1.0 * in2mm;
b.t_fr_nom = 0.44 * in2mm;
b.t_fr_wecsim = 0.44 * in2mm;
b.t_fr_start = 0.44 * in2mm;

% materal thickness of float circumferential gussets (mm)
b.t_fc_min = 0.1 * in2mm;
b.t_fc_max = 1.0 * in2mm;
b.t_fc_nom = 0.44 * in2mm;
b.t_fc_wecsim = 0.44 * in2mm;
b.t_fc_start = 0.44 * in2mm;

% material thickness of float bottom (mm)
b.t_fb_min = 0.1 * in2mm;
b.t_fb_max = 1.0 * in2mm;
b.t_fb_nom = 0.56 * in2mm;
b.t_fb_wecsim = 0.56 * in2mm;
b.t_fb_start = 0.56 * in2mm;

% material thickness of spar radial (mm)
b.t_sr_min = 0.2 * in2mm;
b.t_sr_max = 2.0 * in2mm;
b.t_sr_nom = 1.0 * in2mm;
b.t_sr_wecsim = 1.0 * in2mm;
b.t_sr_start = 1.0 * in2mm;

% material thickness of damping plate support tube radial walls (mm)
b.t_dt_min = 0.2 * in2mm;
b.t_dt_max = 2.0 * in2mm;
b.t_dt_nom = 1.0 * in2mm;
b.t_dt_wecsim = 1.0 * in2mm;
b.t_dt_start = 1.0 * in2mm;

% maximum generator power (kW)
b.P_max_min = 50;
b.P_max_max = 1000;
b.P_max_nom = 286;
b.P_max_wecsim = 286;
b.P_max_start = 286;

% material index (-)
b.M_min = 1;
b.M_max = 3;
b.M_wecsim = 1;
b.M_nom = 1;
b.M_start = 1;

% D_s D_f T_f_2 h_s F_max B_p h_fs_clear]
b.mins_flexible = [false true true true true true true]';
b.maxs_flexible = [true true false false true true true]';
% D_s D_f T_f_2 h_s F_max B_p h_fs_clear thicknesses P_max]
b.mins_flexible = [false true true true true true true true(1,6) true]';
b.maxs_flexible = [true true false false true true true true(1,6) true]';
% if a bound is marked flexible and the bound is active after optimization,
% a warning in gradient_optim will remind you to adjust the bound.

b.X_mins = [b.D_s_min b.D_f_min b.T_f_2_min b.h_s_min b.F_max_min b.B_p_min b.h_fs_clear_min]';
b.X_maxs = [b.D_s_max b.D_f_max b.T_f_2_max b.h_s_max b.F_max_max b.B_p_max b.h_fs_clear_max]';
b.X_mins = [b.D_s_min b.D_f_min b.T_f_2_min b.h_s_min b.F_max_min b.B_p_min b.h_fs_clear_min ...
b.t_ft_min b.t_fr_min b.t_fc_min b.t_fb_min b.t_sr_min b.t_dt_min b.P_max_min]';
b.X_maxs = [b.D_s_max b.D_f_max b.T_f_2_max b.h_s_max b.F_max_max b.B_p_max b.h_fs_clear_max ...
b.t_ft_max b.t_fr_max b.t_fc_max b.t_fb_max b.t_sr_max b.t_dt_max b.P_max_max]';

if strcmpi(mode,'wecsim')
b.X_noms = [b.D_s_wecsim b.D_f_wecsim b.T_f_2_wecsim b.h_s_wecsim b.F_max_wecsim b.B_p_wecsim b.h_fs_clear_wecsim]';
b.X_noms = [b.D_s_wecsim b.D_f_wecsim b.T_f_2_wecsim b.h_s_wecsim b.F_max_wecsim ...
b.B_p_wecsim b.h_fs_clear_wecsim b.t_ft_wecsim b.t_fr_wecsim b.t_fc_wecsim ...
b.t_fb_wecsim b.t_sr_wecsim b.t_dt_wecsim b.P_max_wecsim]';
else
b.X_noms = [b.D_s_nom b.D_f_nom b.T_f_2_nom b.h_s_nom b.F_max_nom b.B_p_nom b.h_fs_clear_nom]';
b.X_noms = [b.D_s_nom b.D_f_nom b.T_f_2_nom b.h_s_nom b.F_max_nom b.B_p_nom b.h_fs_clear_nom ...
b.t_ft_nom b.t_fr_nom b.t_fc_nom b.t_fb_nom b.t_sr_nom b.t_dt_nom b.P_max_nom]';
end
b.X_starts = [b.D_s_start b.D_f_start b.T_f_2_start b.h_s_start b.F_max_start b.B_p_start b.h_fs_clear_start]';

b.X_starts = [b.D_s_start b.D_f_start b.T_f_2_start b.h_s_start b.F_max_start b.B_p_start b.h_fs_clear_start ...
b.t_ft_start b.t_fr_start b.t_fc_start b.t_fb_start b.t_sr_start b.t_dt_start b.P_max_start]';

b.X_start_struct = cell2struct(num2cell(b.X_starts),b.var_names(1:end-1)',1);

Expand All @@ -92,9 +152,11 @@
for i = (i1+1):(i1+14*15)
b.constraint_names{i} = strcat('prevent_slamming',num2str(i-i1));
end
b.constraint_names_pretty = remove_underscores(b.constraint_names);

b.lin_constraint_names = {'spar_natural_freq','float_spar_diam','float_spar_draft',...
'float_spar_tops','float_seafloor','spar_seafloor'};
b.lin_constraint_names_pretty = remove_underscores(b.lin_constraint_names);

[~,idxs_sort] = sort(b.var_names(1:end-1)); % alphabetical design variable indices
idxs_recover = zeros(size(idxs_sort));
Expand Down
43 changes: 29 additions & 14 deletions mdocean/optimization/gradient_optim.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,13 +21,20 @@

% create optimization variables for each of the design variables
sz = [1 1]; % create scalar variables
x1 = optimvar(b.var_names{1}, sz,'LowerBound',b.X_mins(1), 'UpperBound',b.X_maxs(1));
x2 = optimvar(b.var_names{2}, sz,'LowerBound',b.X_mins(2), 'UpperBound',b.X_maxs(2));
x3 = optimvar(b.var_names{3}, sz,'LowerBound',b.X_mins(3), 'UpperBound',b.X_maxs(3));
x4 = optimvar(b.var_names{4}, sz,'LowerBound',b.X_mins(4), 'UpperBound',b.X_maxs(4));
x5 = optimvar(b.var_names{5}, sz,'LowerBound',b.X_mins(5), 'UpperBound',b.X_maxs(5));
x6 = optimvar(b.var_names{6}, sz,'LowerBound',b.X_mins(6), 'UpperBound',b.X_maxs(6));
x7 = optimvar(b.var_names{7}, sz,'LowerBound',b.X_mins(7), 'UpperBound',b.X_maxs(7));
x1 = optimvar(b.var_names{1}, sz,'LowerBound',b.X_mins(1), 'UpperBound',b.X_maxs(1));
x2 = optimvar(b.var_names{2}, sz,'LowerBound',b.X_mins(2), 'UpperBound',b.X_maxs(2));
x3 = optimvar(b.var_names{3}, sz,'LowerBound',b.X_mins(3), 'UpperBound',b.X_maxs(3));
x4 = optimvar(b.var_names{4}, sz,'LowerBound',b.X_mins(4), 'UpperBound',b.X_maxs(4));
x5 = optimvar(b.var_names{5}, sz,'LowerBound',b.X_mins(5), 'UpperBound',b.X_maxs(5));
x6 = optimvar(b.var_names{6}, sz,'LowerBound',b.X_mins(6), 'UpperBound',b.X_maxs(6));
x7 = optimvar(b.var_names{7}, sz,'LowerBound',b.X_mins(7), 'UpperBound',b.X_maxs(7));
x8 = optimvar(b.var_names{8}, sz,'LowerBound',b.X_mins(8), 'UpperBound',b.X_maxs(8));
x9 = optimvar(b.var_names{9}, sz,'LowerBound',b.X_mins(9), 'UpperBound',b.X_maxs(9));
x10 = optimvar(b.var_names{10}, sz,'LowerBound',b.X_mins(10), 'UpperBound',b.X_maxs(10));
x11 = optimvar(b.var_names{11}, sz,'LowerBound',b.X_mins(11), 'UpperBound',b.X_maxs(11));
x12 = optimvar(b.var_names{12}, sz,'LowerBound',b.X_mins(12), 'UpperBound',b.X_maxs(12));
x13 = optimvar(b.var_names{13}, sz,'LowerBound',b.X_mins(13), 'UpperBound',b.X_maxs(13));
x14 = optimvar(b.var_names{14}, sz,'LowerBound',b.X_mins(14), 'UpperBound',b.X_maxs(14));

opts = optimoptions('fmincon', 'Display',display,...
'Algorithm','sqp',...%'interior-point',...
Expand All @@ -42,7 +49,7 @@

% iterate through material choices
for matl = 1%1:2:3 %b.M_min : b.M_max
X = [x1 x2 x3 x4 x5 x6 x7 matl];
X = [x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 matl];
[Xs_opt, objs_opt, flags, probs] = optimize_both_objectives(X,p,b,x0_input,opts,ploton,which_objs);

end
Expand Down Expand Up @@ -111,12 +118,20 @@
probs{i} = problem;

tol = eps(2);
if any(abs(X_opt_raw(b.mins_flexible) - b.X_mins(b.mins_flexible)) < tol) ...
|| any(abs(X_opt_raw(b.maxs_flexible) - b.X_maxs(b.maxs_flexible)) < tol)
warning('Optimization is up against a flexible variable bound, consider changing bounds')
min_active = abs(X_opt_raw - b.X_mins) < tol;
max_active = abs(X_opt_raw - b.X_maxs) < tol;
flexible_min_active = b.mins_flexible & min_active;
flexible_max_active = b.maxs_flexible & max_active;
if any(flexible_min_active) || any(flexible_max_active)
var_opt = b.var_names(1:end-1);
max_string = strcat(var_opt(flexible_max_active)," ");
min_string = strcat(var_opt(flexible_min_active)," ");
msg = ['Optimization is up against a flexible variable bound, consider changing bounds. ',...
'At max: ' horzcat(max_string{:}) ', at min: ' horzcat(min_string{:})];
warning(msg)
end

X_opt = [X_opt_raw; evaluate(X(8),struct())]; % add material back onto design vector
X_opt = [X_opt_raw; evaluate(X(end),struct())]; % add material back onto design vector
[out(1),out(2)] = simulation(X_opt,p); % rerun sim
assert(out(which_obj) == obj_opt) % check correct reordering of X_opt elements

Expand All @@ -131,11 +146,11 @@
end
end
if ploton
table_data = [Xs_opt(1:end-1,:), b.X_mins, b.X_maxs];
table_data = [Xs_opt(1:end-1,:), b.X_mins, b.X_maxs b.X_noms];
objs_opt
flags
array2table(table_data,'RowNames',b.var_names(1:end-1),...
'VariableNames',{'Min LCOE','Min cv','Min bound','Max bound'})
'VariableNames',{'Min LCOE','Min cv','Min bound','Max bound','Nom'})
end

end
14 changes: 10 additions & 4 deletions mdocean/optimization/multiobjective/pareto_curve_heuristics.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,12 @@ function pareto_curve_heuristics()
%[x,fval] = pareto_search();
d=dir("**/pareto_search_results*");
load(d(end).name)

if ~exist('tol','var')
tol = 1e-6;
end
constraint_active_plot(residuals,fval,tol,b)

cols = b.idxs_recover;
X = x(:,cols); % swap indices based on solver generated function
X = [X ones(length(X),1)]; % add extra column for material
Expand Down Expand Up @@ -207,13 +213,13 @@ function pareto_curve_heuristics()

X_pareto_sorted_scaled = X_pareto_sorted ./ repmat(x_best_LCOE,length(idx_sort),1);

X_pareto_sorted_scaled = X_pareto_sorted_scaled(:,1:7); % get rid of material
X_pareto_sorted_scaled = X_pareto_sorted_scaled(:,1:end-1); % get rid of material

windowSize = round(length(idx_sort) * 5/100);
b = (1/windowSize)*ones(1,windowSize);
a = 1;
y = zeros(size(X_pareto_sorted_scaled));
for i=1:7
for i=1:size(X_pareto_sorted_scaled,2)
x = X_pareto_sorted_scaled(:,i);
x_padded = [ones(windowSize,1); x];
yy = filter(b,a,x_padded); % moving average filter
Expand All @@ -233,9 +239,9 @@ function pareto_curve_heuristics()
set(gca,'YMinorGrid','on')

% filtered
cols = {'r:','r--','r-','r-.','b:','b--','b-'};
cols = {'r:','r--','r-','r-.','b:','b--','b-','g:','g--','g-','g-.','g.','g*','b.'};
figure
for i=1:7
for i=1:size(X_pareto_sorted_scaled,2)
semilogy(pct_angle,y(:,i),cols{i})
hold on
end
Expand Down
Loading

0 comments on commit ee986b1

Please sign in to comment.