Merge pull request #29 from symbiotic-engineering/pareto-front-tweak

Pareto front tweak

CITATION.cff

@@ -11,7 +11,7 @@ preferred-citation:
       orcid: "https://orcid.org/0000-0000-0000-0000"
     - family-names: "Haji"
       given-names: "Maha N."
-      orcid: "https://orcid.org/0000-0000-0000-0000"
+      orcid: "https://orcid.org/0000-0002-2953-7253"
   title: "Multidisciplinary Optimization to Reduce Cost and Power Variation of a Wave Energy Converter"
   year: 2022
   conference:

README.md

@@ -11,7 +11,7 @@ in the [2022 ASME IDETC-CIE](https://event.asme.org/IDETC-CIE).
 At this conference, the work was presented at the [DAC-6](https://www.designautomationconference.org/dac-6) session and is publication number 90227.
 The project began as an effort in Cornell course [MAE 5350](https://classes.cornell.edu/browse/roster/FA21/class/MAE/5350).
 
-Citation: `R. McCabe, O. Murphy, and M. N. Haji, “Multidisciplinary Optimization to Reduce Cost and Power Variation of a Wave Energy Converter,” International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, St. Louis, MO, August 14-17, 2022.`
+Citation: R. McCabe, O. Murphy, and M. N. Haji, “Multidisciplinary Optimization to Reduce Cost and Power Variation of a Wave Energy Converter,” *International Design Engineering Technical Conferences & Computers and Information in Engineering Conference*, St. Louis, MO, August 14-17, 2022.
 
 **Authors**
 - Rebecca McCabe, [email protected] (Project lead and point of contact) @rebeccamccabe
@@ -21,7 +21,7 @@ Citation: `R. McCabe, O. Murphy, and M. N. Haji, “Multidisciplinary Optimizati
 **Disclaimer**
 
 The versions of the simulation used in the conference proceedings paper ([v1.2](https://github.com/symbiotic-engineering/MDOcean/releases/tag/v1.2)) and in the conference video ([v1.3](https://github.com/symbiotic-engineering/MDOcean/releases/tag/v1.3)) were later found to have a number of errors. 
-These errors have since been corrected, and the current code is correct to the best of the authors' knowledge, within the limitations of the stated assumptions. 
+These errors have since been corrected, and the current code ([v1.4](https://github.com/symbiotic-engineering/MDOcean/releases/tag/v1.4)) is correct to the best of the authors' knowledge, within the limitations of the stated assumptions. 
 Known areas for improvement are listed as GitHub issues. If you find any additional errors, please let us know.
 
 Errors that have since been fixed include:
@@ -35,6 +35,7 @@ Errors that have since been fixed include:
 - Incorrect scaling of cost with number of WECs
 - Powertrain and transmission loss not accounted for
 - Mass of damping plate support tubes not accounted for
+- Float maximum displacement constraint not being enforced
 
 *Because the results in the official conference proceedings contain these errors, an updated unofficial version 
 of the paper will be posted here in early September 2022, containing corrections.*
@@ -64,4 +65,13 @@ The following packages are used in this code:
 - Statistics and Machine Learning Toolbox
 - Symbolic Math Toolbox
 
-All are required except the symbolic math toolbox, which is used only for code generation and exploratory scripting, not core functionality.
+All are required except the symbolic math toolbox, which is used only for code generation and exploratory scripting, not core functionality.
+
+**Funding Acknowledgement**
+
+This material is based upon work supported by the 
+National Science Foundation Graduate Research Fellowship under 
+Grant No. DGE–2139899, and the Cornell Engineering Fellowship.
+Any opinion, findings, and conclusions or recommendations 
+expressed in this material are those of the authors(s) and do not 
+necessarily reflect the views of the National Science Foundation.

dev/gamma_integral.m

@@ -0,0 +1,138 @@
+syms w t k rho R a T g real positive
+syms x y real
+phi_I = a*w/k * exp(-k*T) * 1i * exp(1i*(w*t - k*x));
+integrand_pressure = rho * diff(phi_I,t);
+integrand_gamma = simplify(integrand_pressure / (a*exp(1i*w*t)));
+
+x_min = -sqrt(R^2-y^2);
+x_max = sqrt(R^2-y^2);
+inner_int = int(integrand_gamma,x,x_min,x_max);
+
+coeffs = rho*w^2*exp(-T*k)/k^2;
+inner_int_without_coeffs = inner_int / coeffs;
+
+y_min = -R;
+y_max = R;
+gamma_without_coeffs = int(inner_int_without_coeffs,y,y_min,y_max);
+gamma_without_coeffs = simplify(rewrite(gamma_without_coeffs,'sin'));
+coeffs = subs(coeffs,k^2,k*w^2/g);
+gamma = abs(gamma_without_coeffs * coeffs)
+taylor(gamma_without_coeffs/k,k,'Order',11)
+
+derivative = diff(gamma_without_coeffs,'k')
+
+%% approximation error
+i = 8:2:14; % last term included
+exc = i+2;  % first term excluded
+w = 1.4;
+k = w^2/9.8;
+R_max = 20;
+order_err = (k*R_max).^exc ./ factorial(exc);
+figure
+plot(i,order_err*100)
+xlabel('exponent of k on last term included')
+ylabel('Approx percent error, 0-100% scale')
+
+%%
+
+
+% syms r theta
+% integrand_gamma_polar = r * subs(integrand_gamma,x,r*cos(theta));
+% inner_int = int(integrand_gamma_polar,r,0,R);
+% gamma = int(inner_int, theta, 0, 2*pi)
+
+% subscript _num is for numerical, so I don't overwrite symbolic vars
+R_num = [3 6 10 12 15 18];
+T_wave = 5:2:13;
+w_num = 2*pi./T_wave;
+[R_mesh,w_mesh] = meshgrid(R_num,w_num);
+k_mesh = w_mesh.^2/9.8;
+approx = 2.88*R_mesh.^2;
+actual = 1./k_mesh .* vpa(subs(gamma_without_coeffs,{R,k},{R_mesh,k_mesh}))
+
+%%
+w_mesh = sqrt(k_mesh*9.8);
+actual=abs(actual);
+error = actual./approx;
+
+figure
+subplot 121
+contourf(R_mesh,k_mesh,approx)
+colorbar
+title('approx')
+hold on
+plot([10 10],[min(k_num) max(k_num)],'w--')
+
+subplot 122
+contourf(R_mesh,k_mesh,actual)
+colorbar
+title('actual')
+hold on
+plot([10 10],[min(k_num) max(k_num)],'w--')
+
+figure
+subplot 121
+contourf(R_mesh,k_mesh,error)
+colorbar
+caxis([0 1])
+title('error ratio')
+hold on
+plot([10 10],[min(k_num) max(k_num)],'w--')
+
+subplot 122
+trial = 1./(k_mesh .* R_mesh); %exp(-w_mesh.^2*0.5
+contourf(R_mesh,k_mesh,trial/2)
+colorbar
+caxis([0 1])
+%%
+% line for each R
+figure
+subplot 121
+for i=1:length(R_num)
+plot(w_num,error(:,i),'DisplayName',num2str(R_num(i)))
+hold on
+end
+title('Error for various radii')
+xlabel('k')
+legend
+
+% line for each k
+subplot 122
+for i=1:length(k_num)
+plot(R_num,error(i,:),'DisplayName',num2str(k_num(i)))
+hold on
+end
+title('Error for various k')
+xlabel('R')
+legend
+
+%%
+% compare numerical integral to actual wamit
+w_wamit = hydro.w;
+gamma = hydro.ex_ma(3,1,:);
+gamma = gamma(:)';
+draft = 2;
+wamit_compare = gamma ./ exp(-w_wamit.^2*draft/9.8)
+mdocean = 314 * exp(w_wamit.^2*draft/9.8);
+
+figure
+plot(w_num,actual(:,R_num==10)-100)
+hold on
+plot(w_num,approx(:,R_num==10))
+plot(w_wamit,mdocean)
+plot(w_wamit,wamit_compare)
+legend('Numerical integral','Approximate Integral (Desmos)','MDOcean original','WAMIT')
+xlim([min(w_num) max(w_num)])
+xlabel('w')
+ylabel('$\frac{\gamma}{\rho g e^{-kT}}$','Interpreter','latex','FontSize',20)
+
+%%
+figure
+plot(w_wamit,gamma)
+hold on
+fudge_factor = 1;
+plot(w_num,actual(:,R_num==10) .* exp(-k_num'*fudge_factor*draft))
+legend('WAMIT','Numerical Integration with 5x fudge')
+xlim([min(w_num) max(w_num)])
+xlabel('w')
+ylabel('\gamma/\rho g')

dev/max_capture_width.m

@@ -0,0 +1,5 @@
+syms T g
+w = 2*pi/T; % definition of angular frequency
+k = w^2/g; % dispersion relation for deep water
+lambda = 2*pi/k; % definition of wavenumber
+CW_max = lambda / (2*pi) % eq 3 https://www.sciencedirect.com/science/article/pii/S0960148115001652#bib50

dev/nominal_c_v.m

@@ -0,0 +1,15 @@
+function c_v = nominal_c_v()
+	filename = 'RM3-CBS.xlsx'; % spreadsheet containing RM3 "actual" power data
+	sheet = 'Performance & Economics'; % name of relevant sheet
+
+	P_matrix = readmatrix(filename,'Range','E97:S110','Sheet',sheet) * 1000;
+	JPD = readmatrix(filename,'Range','E24:S37','Sheet',sheet);
+
+	P_weighted = P_matrix .* JPD / sum(JPD,'all');
+	P_elec = sum(P_weighted(:));
+
+	% coefficient of variance (normalized standard deviation) of power
+	P_var = std(P_matrix(:), JPD(:)) / P_elec;
+	c_v = P_var * 100; % convert to percentage
+
+end

dev/power_matrix_compare.m

@@ -0,0 +1,106 @@
+clear;close all;clc
+
+filename = 'RM3-CBS.xlsx'; % spreadsheet containing RM3 "actual" power data
+
+% inputs
+p = parameters();
+b = var_bounds(p);
+X = [b.X_noms; 1];
+
+% unsaturated power
+actual_mech_unsat = readmatrix(filename,'Range','E73:S86','Sheet','Performance & Economics');
+actual_elec_unsat = actual_mech_unsat * p.pto_eff;
+[~, P_var, ~, ~, ~, ~, ~, ~, P_elec, ~, P_matrix] = simulation(X,p);
+sim_elec_unsat = P_matrix/1000;
+
+% saturated power
+v = validation_inputs();
+p.power_max = v.power_max;
+actual_elec_sat = readmatrix(filename,'Range','E97:S110','Sheet','Performance & Economics');
+[~, P_var, ~, ~, ~, ~, ~, ~, P_elec, ~, P_matrix] = simulation(X,p);
+sim_elec_sat = P_matrix/1000;
+
+% wave resources
+[T,H] = meshgrid(p.T, p.Hs);
+JPD = p.JPD;
+JPD_actual = readmatrix(filename,'Range','E24:S37','Sheet','Performance & Economics');
+wave_resource_raw = 1030 * 9.8^2 / (64*pi) * T .* H.^2 / 1000;
+wave_resource_sheet = readmatrix(filename,'Range','E49:S62','Sheet','Performance & Economics');
+wave_resource_sheet(wave_resource_sheet == 0) = NaN;
+wave_resource_sim = wave_resource_raw;
+wave_resource_sim(JPD == 0) = NaN;
+
+% variable manipulation before plotting
+power_titles = {'Actual: Unsaturated','Simulated: Unsaturated','Actual: Saturated','Simulated: Saturated'};
+pre_JPD_power        = actual_elec_unsat;
+pre_JPD_power(:,:,2) = sim_elec_unsat;
+pre_JPD_power(:,:,3) = actual_elec_sat;
+pre_JPD_power(:,:,4) = sim_elec_sat;
+
+post_JPD_power = pre_JPD_power .* JPD / 100;
+
+diameter = X(1);
+CW = pre_JPD_power ./ wave_resource_raw;    % capture width
+CW_max = 9.8 * T.^2 / (4*pi^2);             % max CW for linear hydrodynamics, axisymmetric body
+CWR = CW / diameter;                        % capture width ratio
+CW_to_CW_max = CW ./ CW_max;
+CW_to_CW_max_zeroed = CW_to_CW_max;
+JPD_rep = repmat(JPD,[1 1 4]);
+CW_to_CW_max_zeroed(JPD_rep == 0) = NaN;
+
+vars = pre_JPD_power;
+vars(:,:,:,2) = post_JPD_power;
+vars(:,:,:,3) = CWR;
+vars(:,:,:,4) = CW_to_CW_max;
+vars(:,:,:,5) = CW_to_CW_max_zeroed;
+
+var_names = {'Unweighted Device Power Matrix (kW)','JPD-Weighted Power Matrix (kW)',...
+    'Capture Width Ratio (-)','Capture Width / Max Capture Width (-)','Capture Width / Max Capture Width (-)'};
+
+% four figures showing weigted and unweighted power and CWR
+for fig = 1:size(vars,4)
+    var = vars(:,:,:,fig);
+    figure
+    for i = 1:4
+        subplot(2,2,i)
+        contourf(T,H,var(:,:,i))
+        xlabel('T (s)')
+        ylabel('Hs (m)')
+        title(power_titles{i})
+        colorbar
+    end
+    sgtitle(var_names{fig})
+
+end
+
+% resource figure
+resource_titles = {'Raw Resource','Zeroed Resource (Actual)','Zeroed Resource (Sim)'};
+wave_resources = wave_resource_raw;
+wave_resources(:,:,2) = wave_resource_sheet;
+wave_resources(:,:,3) = wave_resource_sim;
+figure
+for i = 1:3
+    subplot(1,3,i)
+    contourf(T,H,wave_resources(:,:,i))
+    xlabel('Wave Period T (s)')
+    ylabel('Wave Height Hs (m)')
+    title(resource_titles{i})
+    colorbar
+end
+sgtitle('Wave Resource')
+
+% JPD figure
+JPDs = JPD;
+JPDs(:,:,2) = JPD_actual/100;
+JPD_titles = {'Sim','Actual'};
+figure
+for i = 1:2
+    subplot(1,2,i)
+    contourf(T,H,JPDs(:,:,i))
+    xlabel('Wave Period T (s)')
+    ylabel('Wave Height Hs (m)')
+    colorbar
+    grid on
+    title(JPD_titles{i})
+end
+sgtitle('JPD (-)')

dev/wamit_coeff_plot.m

@@ -0,0 +1,53 @@
+hydro = struct();
+hydro = readWAMIT(hydro,'rm3.out',[]); % function from WECSim
+
+w = hydro.w;
+A = hydro.A(3,3,:);
+A = A(:);
+B = hydro.B(3,3,:);
+B = B(:);
+gamma = hydro.ex_ma(3,1,:);
+gamma = gamma(:);
+
+r = 10;
+fudge_factor =  5;
+draft = 2;
+g = 9.8;
+rho_w = 1000;
+
+k = w.^2 / g;
+V_g = g ./(2*w);
+
+A_MDOcean = ones(size(w))* 1/2 * 4/3 * pi * r^3 * 0.63;
+r_k_term = r^2 - 1/8 * k.^2 * r^4 + 1/192 * k.^4 * r^6 - 1/9216 * k.^6 * r^8;
+gamma_MDOcean   = pi * exp(-k * draft * fudge_factor) .* r_k_term; 
+B_MDOcean = k/2 .* gamma_MDOcean.^2;
+%% coeff comparison validation figure
+figure
+plot(w,A,'--',w,B,'--',w,gamma,'--','LineWidth',3,'HandleVisibility','off')
+hold on
+set(gca,"ColorOrderIndex",1)
+plot(w,A_MDOcean,w,B_MDOcean,w,gamma_MDOcean)
+xlim([0 1.5])
+ylim([0 2500])
+plot(0,0,'k-',0,0,'k--') % dummy plot so I can get 2 extra legend entries
+legend('Added Mass A/\rho','Radiation Damping B/(\rho\omega)','Excitation Force \gamma/(\rhog)','Simulation (Analytical)','Actual (WAMIT BEM)')
+title('Normalized Hydrodynamic Coefficients')
+xlabel('Wave frequency \omega (rad/s)')
+improvePlot
+
+%% R^2 table
+
+R_A = corrcoef(A, A_MDOcean);
+R_B = corrcoef(B, B_MDOcean);
+R_G = corrcoef(gamma, gamma_MDOcean);
+
+R2_A = R_A(1,2)^2
+R2_B = R_B(1,2)^2
+R2_G = R_G(1,2)^2
+
+%% check B formula
+figure
+plot(w,B./gamma.^2*10000, w,w.^2/(2*9.8)*10000,'--')
+legend('B/\gamma^2*\rho g^2/\omega', '\omega^2/2g')
+xlim([0 1.5])