diff --git a/.gitignore b/.gitignore index 64b8c976..1770b5b5 100644 --- a/.gitignore +++ b/.gitignore @@ -45,6 +45,8 @@ test-results/ # figures *.fig *.png +# except intentional figures +!**/non_matlab_figs/*.png # mat data files *.mat diff --git a/dev/dynamics/multibody/multibody_fudge_factor.m b/dev/dynamics/multibody/multibody_fudge_factor.m new file mode 100644 index 00000000..caec4692 --- /dev/null +++ b/dev/dynamics/multibody/multibody_fudge_factor.m @@ -0,0 +1,22 @@ +p = parameters(); +b = var_bounds(); +X = [b.X_noms; 1]; +wecsim_filename = 'wecsim_sparcd0_floatcd0_14e825b_01a25103-1490-4189-b376-76cbdf616634'; +report = true; +override = false; +power_matrix_compare(X, p, wecsim_filename, report, override) + +% add breakpoint inside power_matrix_compare and do the following in +% command window: (see 1/27/25 slides for screenshots) +sp = sim.power_mech_unsat; +ap = actual.power_mech_unsat; +ratio = ap ./ sp + +figure; plot(T,ratio) +ylabel('RM3 report power / wecsim singlebody power') +xlabel('T') + +y=ratio(H==5.75,:); +figure; plot(T,y) + +% then do apps > curve fitter with rational function: num degree 0, den degree 2 \ No newline at end of file diff --git a/dev/econ/design_cost_scaling.m b/dev/econ/design_cost_scaling.m new file mode 100644 index 00000000..2463da42 --- /dev/null +++ b/dev/econ/design_cost_scaling.m @@ -0,0 +1,50 @@ +%%%% PTO cost %%%%%% +% PTO cost scale with N_WEC (from curve fit) +% C_PTO_per_wec = C_PTO_const + C_PTO_coeff * N_WEC ^ -C_PTO_alpha +C_PTO_const = 0.280e6; +C_PTO_coeff = 0.344e6; +C_PTO_alpha = 0.206; + +% cost that scales with force: generator + bearing + mounting +cost_force_N1 = 25740+88411+28600; +cost_force_N100 = 20037+66308+5720; + +F_max_nom = 0.926e6; +C_force_coeff = (cost_force_N1 - cost_force_N100) / (F_max_nom * (1 - 100^(-C_PTO_alpha))) +C_force_const = (cost_force_N1 / F_max_nom) - C_force_coeff + +% cost that scales with neither force nor power: riser cable and controller +cost_neither_N1 = 88000+5644; +cost_neither_N100 = 88000+5000; + +C_neither_coeff = (cost_neither_N1 - cost_neither_N100) / ((1 - 100^(-C_PTO_alpha))) +C_neither_const = (cost_neither_N1) - C_neither_coeff + +P_max_nom = 286000; +C_power_coeff = (C_PTO_coeff - C_neither_coeff - C_force_coeff*F_max_nom) / P_max_nom +C_power_const = (C_PTO_const - C_neither_const - C_force_const*F_max_nom) / P_max_nom + + +cost_force_pred = (C_force_const + C_force_coeff*[1 100].^-C_PTO_alpha) * F_max_nom; +cost_force_act = [cost_force_N1 cost_force_N100]; + +cost_neither_pred = C_neither_const + C_neither_coeff*[1 100].^-C_PTO_alpha; +cost_neither_act = [cost_neither_N1 cost_neither_N100]; + +cost_pto_pred = C_PTO_const + C_PTO_coeff*[1 100].^-C_PTO_alpha +cost_power_pred = cost_pto_pred - cost_neither_pred - cost_neither_pred + +%%%%%% non design cost %%%%%%%% +clc +N_WEC = [1 10 50 100]; +development = [4553389 8773812 11003159 10820060] ./ N_WEC +infrastructure = [990000 4860000 7566000 17310000] ./ N_WEC +mooring = [524775 4722975 23614875 47229750] ./ N_WEC +profitmargin = [356252 2561152 11323295 21921723] ./ N_WEC +installation = [5908552 9081973 21531225 37859591] ./ N_WEC +decommissioning = 0;%installation; % decomissioning cost not used in LCOE +contingency = [1589545 5561144 18827150 35435836] ./ N_WEC + +non_design_capex = development + infrastructure + mooring + profitmargin + installation + decommissioning + contingency; + + diff --git a/dev/lcoe_trend.m b/dev/econ/lcoe_trend.m similarity index 100% rename from dev/lcoe_trend.m rename to dev/econ/lcoe_trend.m diff --git a/dev/force_sat/quadratic.m b/dev/force_sat/quadratic.m index 9c784d8a..ae49b368 100644 --- a/dev/force_sat/quadratic.m +++ b/dev/force_sat/quadratic.m @@ -68,7 +68,8 @@ matrix = simplify(matrix,'IgnoreAnalyticConstraints',true) latex(matrix) -bases_rk = bases .* [r_k^2, r_k^2*r_b, r_k^2*r_b^2, r_k*r_b^2, r_b^2]; +rkb = [r_k^2, r_k^2*r_b, r_k^2*r_b^2, r_k*r_b^2, r_b^2]; +bases_rk = bases .* rkb; pretty(bases_rk') latex(bases_rk) @@ -77,6 +78,51 @@ soln_f = solve(eqn_m_1,f); assert(soln_f == 1) +%% replace rk with other variables +% notebook p97 2-1-25 +syms omega_u_star zeta_u star real positive + +rk_new = star^2/(star^2-1); +rb_new = 1/(1-zeta_u/(zeta*star)); +w_rat_new = omega_u_star / star; + +rkb_new = subs(rkb,[r_k,r_b],[rk_new, rb_new]); +pretty(rkb_new.') +bases_new = subs(bases,w_ratio,w_rat_new); +bases_rk_new = bases_new .* rkb_new; +scale = (star^2-1)^2 * (star*zeta)^-2 * (zeta_u - star*zeta)^2; +pretty(simplify(bases_rk_new.' * scale,'IgnoreAnalyticConstraints',true)) + +%% check that parameterization used in paper equals "new" just above +% notebook p99 2-1-25 +matrix_check = [-4 4/f^2 1/f^2 -1 0 0; + 8 0 0 2 8 2; + -4 -4 -1 -1 -8 -2]; +bases_unsquared_check = [omega_u_star * (zeta_u - star*zeta),... + omega_u_star*star*zeta,... + omega_u_star^2 - star^2,... + star^2-1]; +bases_squared_check = [bases_unsquared_check(1)^2;... + bases_unsquared_check(2)^2;... + bases_unsquared_check(3)^2;... + bases_unsquared_check(4)^2;... + bases_unsquared_check(1)*bases_unsquared_check(2);... + bases_unsquared_check(3)*bases_unsquared_check(4)]; + +bases_should_match = [bases_unsquared_check(1)^2;... % a^2 + -bases_unsquared_check(1) * bases_unsquared_check(2);... % -ab + bases_unsquared_check(3)^2 + 4*bases_unsquared_check(2)^2;... % c^2+4b^2 + bases_unsquared_check(3) * bases_unsquared_check(4);... % cd + bases_unsquared_check(4)^2]; % d^2 + +should_be_1 = simplify( bases_should_match ./ (bases_rk_new.' * scale) ) +assert(all(should_be_1 == 1)) + +%% and check that it equals the original (r_k, r_b, quad) +quad_new = subs(quad,[r_k,r_b,w_ratio],[rk_new,rb_new,w_rat_new]); +should_be_1_also = simplify( quad_new.' ./ (matrix_check*bases_squared_check/scale) ) +assert(all(should_be_1_also == 1)) + %% optimal powertrain case bases_rk = subs(bases_rk,{w_ratio,r_b},{1,2}).'; diff --git a/mdocean/inputs/parameters.m b/mdocean/inputs/parameters.m index dccac0da..fa23a6ee 100644 --- a/mdocean/inputs/parameters.m +++ b/mdocean/inputs/parameters.m @@ -42,7 +42,7 @@ table("rho_w","\rho_w",{1000},"site",false,"water density (kg/m3)"); ...%table("mu","\mu",1e-3,"site",false,"dynamic viscosity of water (Pa s)"); table("g","g",{g},"site",false,"acceleration of gravity (m/s2)"); - table("h","h",{100},"site",true,"water depth (m)"); + table("h","h",{45},"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,"significant wave height (m)"); @@ -61,36 +61,41 @@ table("rho_m","\rho_m",{[7850 2400 7900]},"structures",true,"material density (kg/m3)"); table("E","E",{[200e9, 5000*sqrt(4.5*ksi2pa), 200e9]},"structures",true,... "young's modulus (Pa)"); - table("cost_m","cost_m",{[4.28, 125/yd2m^3/2400, 1.84/lb2kg]},"economics",true,"material cost ($/kg)"); + table("cost_perkg_mult","cost_m",{[4.28, 125/yd2m^3/2400, 1.84/lb2kg]/4.28},"economics",true,"material cost ($/kg)"); ...% RM3 CBS sheet 1.4 average of cells F21, F34, F46, F56 ...% https://www.concretenetwork.com/concrete-prices.html ...% https://agmetalminer.com/metal-prices/ table("nu","\nu",{[0.36 0 0.29]},"structures",true,"Poisson's ratio (-)"); + table("FOS_min","FOS_{min}",{1.5},"structures",true,"minimum FOS (-)"); + ... + ...% Thicknesses and structures: float + table("t_f_t_over_t_f_b","t_{f,t}/t_{f,b}",{0.50/0.56},"structures",true,... + "float top to bottom thickness ratio (-)"); + table("t_f_r_over_t_f_b","t_{f,r}/t_{f,b}",{0.44/0.56},"structures",true,... + "float radial to bottom thickness ratio (-)"); + table("t_f_c_over_t_f_b","t_{f,c}/t_{f,b}",{0.44/0.56},"structures",true,... + "float circumferential to bottom thickness ratio (-)"); + table("D_f_tu","D_{f,tu}",{20 * in2m},"structures",true,... + "float support tube diameter (m)"); % 24 in p156 report, 20 in cad + table("t_f_tu","t_{f,tu}",{.5 * in2m},"structures",true,... + "float support tube thickness (m)"); + table("w_over_h_stiff_f","w_{stiff,f}/h_{stiff,f}",{1/16},"structures",... + true,"float stiffener width to height ratio (-)"); + table("num_sections_f","N_{sect}",{12},"structures",false,"number of float sections (-)"); ... - ...% 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_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)"); + ... Thicknesses and structures: damping plate + 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,... "angle from horizontal of damping plate support tubes (rad)"); - table("FOS_min","FOS_{min}",{1.5},"structures",true,"minimum FOS (-)"); - table("D_f_tu","D_{f,tu}",{20 * in2m},"structures",true,"float support tube diameter (m)"); % 24 in p156 report, 20 in cad - table("t_f_tu","t_{f,tu}",{.5 * in2m},"structures",true,"float support tube thickness (m)"); - table("h_stiff_f","h_{stiff,f}",{16 * in2m},"structures",true,"float stiffener height (m)"); - table("w_stiff_f","w_{stiff,f}",{1 * in2m},"structures",true,"float stiffener width (m)"); - table("h_stiff_d","h_{stiff,d}",{[12.5 .5 22 1]*in2m},"structures",true,"damping plate stiffener height (m)"); - table("w_stiff_d","w_{stiff,d}",{[.5 10 1 12] * in2m},"structures",true,"damping plate stiffener width (m)"); - table("num_sections","N_{sect}",{12},"structures",false,"number of float sections (-)"); - table("FOS_mult_d","FOS_{mult,d}",{7.5},"structures",true,"damping plate factor of safety multiplier (-)"); + table("h_over_h1_stiff_d","h_{stiff,d}/h_{1,stiff,d}",{[12.5 0.5 22 1]/22},... + "structures",true,"damping plate stiffener height (m)"); + table("w_over_h1_stiff_d","w_{stiff,d}/h_{1,stiff,d}",{[.5 10 1 12]/22},... + "structures",true,"damping plate stiffener width (m)"); + table("FOS_mult_d","FOS_{mult,d}",{7.5},"structures",true,... + "damping plate factor of safety multiplier (-)"); table("num_terms_plate","N_{plate}",{100},"structures",false,... "number of terms for damping plate concentrated load (-)"); table("radial_mesh_plate","N_{r,plate}",{20},"structures",false,... @@ -99,22 +104,24 @@ "number of damping plate stiffeners (-)"); ... ...% Economics - table("m_scale","m_{scale}",{1.1},"economics",false,... - "factor to account for mass of neglected stiffeners (-)"); + table("m_scale","m_{scale}",{1.10},"economics",false,... + "factor to account for mass of neglected stiffeners and device access (-)"); + table("power_scale","P_{scale}",{85.9/117.7},"economics",false,... + "factor to scale power for validation tuning (-)") table("FCR","FCR",{0.113},"economics",true,... "fixed charge rate (-), see RM3 report p63"); table("N_WEC","N_{WEC}",{100},"economics",true,"number of WECs in array (-)"); - table("LCOE_max","LCOE_{max}",{.5},"economics",true,"maximum LCOE ($/kWh)"); - table("avg_power_min","P_{avg,elec,min}",{100},"economics",true,... - "minimum average electrical power (W)"); % set to a negative number (not zero) to disable constraint + table("LCOE_max","LCOE_{max}",{.75},"economics",true,"maximum LCOE ($/kWh)"); +...% table("avg_power_min","P_{avg,elec,min}",{100},"economics",true,... +...% "minimum average electrical power (W)"); % set to a negative number (not zero) to disable constraint table("eff_array","\eta_{array}",{0.95*0.98},"economics",true,... "array availability and transmission efficiency (-)"); - table("cost_perN", "$/N",{0.1656},"economics",true,... - "cost per Newton ($/N), from https://doi.org/10.1016/j.ifacol.2022.10.531"); - table("cost_perW", "$/W",{1},"economics",true,... - "cost per Watt ($/W)"); + table("cost_perN_mult", "$/N",{1},"economics",true,... + "cost per Newton multiplier (-), 0.1656 from https://doi.org/10.1016/j.ifacol.2022.10.531"); + table("cost_perW_mult", "$/W",{1},"economics",true,... + "cost per Watt multiplier (-)"); ... - ...% Geometric ratios + ...% Geometric ratios of bulk dimensions table("D_d_min","D_{d,min}",{30},"geometry",true,... "minimum damping plate diameter"); table("D_d_over_D_s","D_d/D_s",{30/6},"geometry",true,... @@ -133,14 +140,13 @@ true,"ratio of float inner diameter to spar diameter (-)") ... ...% Dynamics: device parameters - table("C_d_float","C_{d,float}",{0},"dynamics",true,"coefficient of drag for float"); - table("C_d_spar","C_{d,spar}",{0},"dynamics",true,"spar coefficient of drag"); - table("power_max","power_{max}",{Inf},"dynamics",true,"maximum power (W)"); + table("C_d_float","C_{d,float}",{1},"dynamics",true,"coefficient of drag for float"); + table("C_d_spar","C_{d,spar}",{1},"dynamics",true,"spar coefficient of drag"); table("eff_pto","\eta_{pto}",{0.80},"dynamics",true,"PTO efficiency (-)"); ... ...% Dynamics: simulation type table("control_type","control type",{'damping'},"dynamics",false,... - "reactive or constant impedance or damping"); + "reactive or damping"); table("use_MEEM","use_MEEM",{true},"dynamics",false,... "whether to use MEEM for hydro coeffs (boolean)"); table("use_multibody","use_multibody",{false},"dynamics",false,... diff --git a/mdocean/inputs/validation/WEC-Sim/readWAMIT.m b/mdocean/inputs/validation/WEC-Sim/readWAMIT.m index 84baa908..e77234f4 100644 --- a/mdocean/inputs/validation/WEC-Sim/readWAMIT.m +++ b/mdocean/inputs/validation/WEC-Sim/readWAMIT.m @@ -34,7 +34,7 @@ F = b+1; end -p = waitbar(0,'Reading WAMIT output file...'); % Progress bar +%p = waitbar(0,'Reading WAMIT output file...'); % Progress bar e = 0; if isempty(exCoeff)==1; exCoeff = 'diffraction'; end % 'diffraction' or 'haskind' @@ -54,8 +54,8 @@ hydro(F).Nb = 0; % Number of bodies k = 0; for n = 1:N - if ((isempty(strfind(raw{n},'Input from Geometric Data File'))==0 | isempty(strfind(raw{n},'N='))==0) &... - (isempty(strfind(raw{n},'.GDF'))==0 | isempty(strfind(raw{n},'.gdf'))==0)) + if ((isempty(strfind(raw{n},'Input from Geometric Data File'))==0 || isempty(strfind(raw{n},'N='))==0) &&... + (isempty(strfind(raw{n},'.GDF'))==0 || isempty(strfind(raw{n},'.gdf'))==0)) k = k+1; tmp = strsplit(raw{n},{' ','.'}); tmp(cellfun('isempty',tmp)) = []; @@ -105,8 +105,8 @@ hydro(F).w(hydro(F).Nf) = 2*pi/tmp{1}; % Wave frequencies end i = 4; - while (isempty(strfind(raw{n+i},'*******************************'))~=0 & ... - isempty(strfind(raw{n+i},'EXCITING FORCES AND MOMENTS'))~=0 & ... + while (isempty(strfind(raw{n+i},'*******************************'))~=0 && ... + isempty(strfind(raw{n+i},'EXCITING FORCES AND MOMENTS'))~=0 && ... isempty(strfind(raw{n+i},'RESPONSE AMPLITUDE OPERATORS'))~=0) tmp = textscan(raw{n+i},'%f'); if T==0 @@ -122,11 +122,11 @@ i = i+1; end end - if ((isempty(strfind(raw{n},'HASKIND EXCITING FORCES AND MOMENTS'))==0 & ... - strcmp(exCoeff,'haskind')==1) |... - (isempty(strfind(raw{n},'DIFFRACTION EXCITING FORCES AND MOMENTS'))==0 & ... - strcmp(exCoeff,'diffraction')==1) |... - (isempty(strfind(raw{n},'RESPONSE AMPLITUDE OPERATORS'))==0 & ... + if ((isempty(strfind(raw{n},'HASKIND EXCITING FORCES AND MOMENTS'))==0 && ... + strcmp(exCoeff,'haskind')==1) ||... + (isempty(strfind(raw{n},'DIFFRACTION EXCITING FORCES AND MOMENTS'))==0 && ... + strcmp(exCoeff,'diffraction')==1) ||... + (isempty(strfind(raw{n},'RESPONSE AMPLITUDE OPERATORS'))==0 && ... strcmp(exCoeff,'rao')==1)) hydro(F).Nh = 0; % Number of wave headings i = n+1; @@ -135,11 +135,11 @@ tmp = textscan(raw{i}(find(raw{i}==':')+1:end),'%f'); hydro(F).theta(hydro(F).Nh) = tmp{1}; % Wave headings i = i+2; - while (isempty(strfind(raw{i},'*******************************'))~=0 & ... - isempty(strfind(raw{i},'EXCITING FORCES AND MOMENTS'))~=0 & ... - isempty(strfind(raw{i},'RESPONSE AMPLITUDE OPERATORS'))~=0 & ... - isempty(strfind(raw{i},'SURGE, SWAY & YAW DRIFT FORCES (Momentum Conservation)'))~=0 & ... - isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Control Surface)'))~=0 & ... + while (isempty(strfind(raw{i},'*******************************'))~=0 && ... + isempty(strfind(raw{i},'EXCITING FORCES AND MOMENTS'))~=0 && ... + isempty(strfind(raw{i},'RESPONSE AMPLITUDE OPERATORS'))~=0 && ... + isempty(strfind(raw{i},'SURGE, SWAY & YAW DRIFT FORCES (Momentum Conservation)'))~=0 && ... + isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Control Surface)'))~=0 && ... isempty(strfind(raw{i},'Wave Heading'))~=0) tmp = textscan(raw{i},'%f'); ma = tmp{1}(2); % Magnitude of exciting force @@ -166,17 +166,17 @@ tmp = textscan(raw{i}(find(raw{i}==':')+1:end),'%f'); hydro(F).theta(hydro(F).Nh) = tmp{1}(1); % Wave headings i = i+2; - while (isempty(strfind(raw{i},'*******************************'))~=0 & ... - isempty(strfind(raw{i},'EXCITING FORCES AND MOMENTS'))~=0 & ... - isempty(strfind(raw{i},'RESPONSE AMPLITUDE OPERATORS'))~=0 & ... - isempty(strfind(raw{i},'SURGE, SWAY & YAW DRIFT FORCES (Momentum Conservation)'))~=0 & ... - isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Control Surface)'))~=0 & ... - isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Pressure Integration)'))~=0 & ... + while (isempty(strfind(raw{i},'*******************************'))~=0 && ... + isempty(strfind(raw{i},'EXCITING FORCES AND MOMENTS'))~=0 && ... + isempty(strfind(raw{i},'RESPONSE AMPLITUDE OPERATORS'))~=0 && ... + isempty(strfind(raw{i},'SURGE, SWAY & YAW DRIFT FORCES (Momentum Conservation)'))~=0 && ... + isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Control Surface)'))~=0 && ... + isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Pressure Integration)'))~=0 && ... isempty(strfind(raw{i},'Wave Heading'))~=0) tmp = textscan(raw{i},'%f'); ma = tmp{1}(2); % Magnitude of exciting force ph = deg2rad(tmp{1}(3)); % Phase of exciting force - if isnan(ma) ==1 + if isnan(ma) ma = 0; ph = 0; end @@ -196,17 +196,17 @@ tmp = textscan(raw{i}(find(raw{i}==':')+1:end),'%f'); hydro(F).theta(hydro(F).Nh) = tmp{1}(1); % Wave headings i = i+2; - while (isempty(strfind(raw{i},'*******************************'))~=0 & ... - isempty(strfind(raw{i},'EXCITING FORCES AND MOMENTS'))~=0 & ... - isempty(strfind(raw{i},'RESPONSE AMPLITUDE OPERATORS'))~=0 & ... - isempty(strfind(raw{i},'SURGE, SWAY & YAW DRIFT FORCES (Momentum Conservation)'))~=0 & ... - isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Control Surface)'))~=0 & ... - isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Pressure Integration)'))~=0 & ... + while (isempty(strfind(raw{i},'*******************************'))~=0 && ... + isempty(strfind(raw{i},'EXCITING FORCES AND MOMENTS'))~=0 && ... + isempty(strfind(raw{i},'RESPONSE AMPLITUDE OPERATORS'))~=0 && ... + isempty(strfind(raw{i},'SURGE, SWAY & YAW DRIFT FORCES (Momentum Conservation)'))~=0 && ... + isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Control Surface)'))~=0 && ... + isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Pressure Integration)'))~=0 && ... isempty(strfind(raw{i},'Wave Heading'))~=0) tmp = textscan(raw{i},'%f'); ma = tmp{1}(2); % Magnitude of exciting force ph = deg2rad(tmp{1}(3)); % Phase of exciting force - if isnan(ma) ==1 + if isnan(ma) ma = 0; ph = 0; end @@ -226,17 +226,17 @@ tmp = textscan(raw{i}(find(raw{i}==':')+1:end),'%f'); hydro(F).theta(hydro(F).Nh) = tmp{1}(1); % Wave headings i = i+2; - while (isempty(strfind(raw{i},'*******************************'))~=0 & ... - isempty(strfind(raw{i},'EXCITING FORCES AND MOMENTS'))~=0 & ... - isempty(strfind(raw{i},'RESPONSE AMPLITUDE OPERATORS'))~=0 & ... - isempty(strfind(raw{i},'SURGE, SWAY & YAW DRIFT FORCES (Momentum Conservation)'))~=0 & ... - isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Control Surface)'))~=0 & ... - isempty(strfind(raw{i},'SURGE, SWAY & YAW DRIFT FORCES (Momentum Conservation)'))~=0 & ... + while (isempty(strfind(raw{i},'*******************************'))~=0 && ... + isempty(strfind(raw{i},'EXCITING FORCES AND MOMENTS'))~=0 && ... + isempty(strfind(raw{i},'RESPONSE AMPLITUDE OPERATORS'))~=0 && ... + isempty(strfind(raw{i},'SURGE, SWAY & YAW DRIFT FORCES (Momentum Conservation)'))~=0 && ... + isempty(strfind(raw{i},'SURGE, SWAY, HEAVE, ROLL, PITCH & YAW DRIFT FORCES (Control Surface)'))~=0 && ... + isempty(strfind(raw{i},'SURGE, SWAY & YAW DRIFT FORCES (Momentum Conservation)'))~=0 && ... isempty(strfind(raw{i},'Wave Heading'))~=0) tmp = textscan(raw{i},'%f'); ma = tmp{1}(2); % Magnitude of exciting force ph = deg2rad(tmp{1}(3)); % Phase of exciting force - if isnan(ma) ==1 + if isnan(ma) ma = 0; ph = 0; end @@ -249,7 +249,10 @@ end d = floor(10*n/N); % Update progress bar every 10%, otherwise slows computation - if d>e waitbar(n/N); e = d; end + if d>e + %waitbar(n/N); + e = d; + end end %% Scattering Force @@ -358,5 +361,5 @@ %% hydro = normalizeBEM(hydro); % For WAMIT this just sorts the data, if neccessary -close(p); +%close(p); end \ No newline at end of file diff --git a/mdocean/inputs/validation/validation_inputs.m b/mdocean/inputs/validation/validation_inputs.m index 3e47ab56..a02e7acf 100644 --- a/mdocean/inputs/validation/validation_inputs.m +++ b/mdocean/inputs/validation/validation_inputs.m @@ -14,6 +14,10 @@ 'LCOE', [4.48, 1.45, 0.85, 0.76], ... ... % B98:F98 in CBS Report Tables ... % 4.48 is estimated from Fig 5-33 p175 RM3 report + 'capex_struct', 1e6*[2.94 20.67 91.55 177.93] ./ [1 10 50 100], ... + ... % CBS (Total), row 24 + 'capex_PTO', 1e6*[0.62 4.94 21.68 41.28] ./ [1 10 50 100],... + ... % CBS (Total), row 29 'power_avg', 85.9e3, ... % S14 in CBS Performance & Economics 'power_max', 286e3, ... % S15 in CBS Performance & Economics 'force_heave', 8500e3,... % p156 RM3 report @@ -22,12 +26,17 @@ ... % instead of D when calculating I 'c_v', nominal_c_v()); % calculated from data in CBS Performance & Economics +RM3_report.capex_design = RM3_report.capex_struct + RM3_report.capex_PTO; +RM3_report.J_capex_design = RM3_report.capex_design / 1e6; + RM3_wecsim = struct(... 'vol_f', 725.833, ...% submerged volume from WAMIT RM3.out 'vol_s', 886.687, ...% submerged volume from WAMIT RM3.out 'CB_f', 1.2929, ... % center of buoyancy below waterline (MDOcean F24 debugging slide 25) 'CG_f', 0.2833, ... % center of gravity below waterline (MDOcean F24 debugging slide 25) + 'power_avg', NaN, ... % fixme input power from wecsim power matrix 'LCOE', NaN(1,4), ... % fixme input LCOE here from wecsim power matrix and MDOcean econ + 'J_capex_design', NaN(1,4), ... % fixme this shouldn't be needed here since it just comes from MDOcean 'c_v', NaN ); % fixme input c_v here from wecsim power matrix if strcmpi(mode,'wecsim') diff --git a/mdocean/inputs/var_bounds.m b/mdocean/inputs/var_bounds.m index bf3d2661..05f62326 100644 --- a/mdocean/inputs/var_bounds.m +++ b/mdocean/inputs/var_bounds.m @@ -6,12 +6,17 @@ mode = ''; end -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) +% design variables +b.var_names = {'D_s','D_f','T_f_2','h_s','h_fs_clear','F_max','P_max',... + 't_fb','t_sr','t_d','h_stiff_f','h1_stiff_d','M'}; +b.var_names_pretty = {'D_s','D_f','T_{f,2}','h_s','h_{fs,clear}','F_{max}','P_{max}',... + 't_{fb}','t_{sr}','t_d','h_{stiff,f}','h_{1,stiff,d}','M'}; +b.var_descs = {'Spar diameter','Float diameter','Float draught','Spar height',... + 'Float-spar clearance height','Maximum force','Maximum power',... + 'Float bottom thickness','Spar radial thickness','Damping plate thickness',... + 'Float stiffener height','Damping plate stiffener height','Material index'}; + +% diameter of spar (m) b.D_s_min = 0; b.D_s_max = 30; b.D_s_nom = 6; @@ -39,20 +44,6 @@ b.h_s_wecsim = 37.9; b.h_s_start = 37.9; -% maximum powertrain force (MN) -b.F_max_min = 0.01; -b.F_max_max = 100; -b.F_max_wecsim = 100; -b.F_max_nom = 100; -b.F_max_start = 5; - -% powertrain damping (MN / (m/s)) -b.B_p_min = .1; -b.B_p_max = 50; -b.B_p_wecsim = 10; -b.B_p_nom = 10; -b.B_p_start = 0.5; - % vertical clearance between float tube support and spar (m) b.h_fs_clear_min = .01; b.h_fs_clear_max = 10; @@ -60,30 +51,23 @@ 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; +% maximum powertrain force (MN) +b.F_max_min = 0.01; +b.F_max_max = 100; +b.F_max_wecsim = 100; +b.F_max_nom = 1; +b.F_max_start = 5; -% 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; +% maximum generator power (kW) +b.P_max_min = 1; +b.P_max_max = 2000; +b.P_max_nom = 286; +b.P_max_wecsim = 286; +b.P_max_start = 286; +in2mm = 25.4; % material thickness of float bottom (mm) -b.t_fb_min = 0.1 * in2mm; +b.t_fb_min = 0.05 * in2mm; b.t_fb_max = 1.0 * in2mm; b.t_fb_nom = 0.56 * in2mm; b.t_fb_wecsim = 0.56 * in2mm; @@ -96,19 +80,27 @@ 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 thickness of damping plate (mm) +b.t_d_min = 0.05 * in2mm; +b.t_d_max = 2.0 * in2mm; +b.t_d_nom = 1.0 * in2mm; +b.t_d_wecsim = 1.0 * in2mm; +b.t_d_start = 1.0 * in2mm; + +in2m = in2mm / 1000; +% float stiffener height (m) +b.h_stiff_f_min = 0; +b.h_stiff_f_max = 2; +b.h_stiff_f_nom = 16 * in2m; +b.h_stiff_f_wecsim = 16 * in2m; +b.h_stiff_f_start = 16 * in2m; + +% damping plate stiffener height +b.h1_stiff_d_min = 0; +b.h1_stiff_d_max = 2; +b.h1_stiff_d_nom = 22 * in2m; +b.h1_stiff_d_wecsim = 22 * in2m; +b.h1_stiff_d_start = 22 * in2m; % material index (-) b.M_min = 1; @@ -117,40 +109,46 @@ b.M_nom = 1; b.M_start = 1; - % 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]'; + % D_s D_f T_f_2 h_s h_fs_clear F_max P_max thicknesses stiffener-heights ] +b.mins_flexible = [false true true true false true true true(1,3) false(1,2)]'; +b.maxs_flexible = [true true false false true true true true(1,3) true(1,2) ]'; % 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.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]'; - +% create vectors of mins, maxs, noms, and starts +n_dv = length(b.var_names)-1; +[X_mins,X_maxs,X_noms,X_noms_wecsim,X_starts] = deal(zeros(n_dv,1)); +for i=1:n_dv + dv_name = b.var_names{i}; + X_mins(i) = b.([dv_name '_min']); + X_maxs(i) = b.([dv_name '_max']); + X_noms(i) = b.([dv_name '_nom']); + X_noms_wecsim(i) = b.([dv_name '_wecsim']); + X_starts(i) = b.([dv_name '_start']); +end +b.X_mins = X_mins; +b.X_maxs = X_maxs; 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.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]'; + b.X_noms = X_noms_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.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]'; + b.X_noms = X_noms; 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.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_starts = X_starts; b.X_start_struct = cell2struct(num2cell(b.X_starts),b.var_names(1:end-1)',1); +% constraints b.constraint_names = {'float_too_heavy','float_too_light','spar_too_heavy','spar_too_light',... 'stability','FOS_float_max','FOS_float_fatigue',... - 'float_top_ratio','float_radial_ratio','float_circ_ratio',... % placeholders for now 'FOS_col_max','FOS_col_fatigue','FOS_plate_max','FOS_plate_fatigue',... 'FOS_col_local_max','FOS_col_local_fatigue',... 'pos_power','LCOE_max','irrelevant_max_force',... - 'spar_height_up','spar_height_down','linear_theory'}; + 'spar_height_up','spar_height_down','float_spar_hit',... + 'linear_theory'}; i1 = length(b.constraint_names); -for i = (i1+1):(i1+14*15) +JPD_size = 14*15; +storm_size = 7; +for i = (i1+1):(i1+JPD_size+storm_size) b.constraint_names{i} = strcat('prevent_slamming',num2str(i-i1)); end b.constraint_names_pretty = remove_underscores(b.constraint_names); @@ -159,18 +157,23 @@ 'float_spar_tops','float_seafloor','spar_seafloor'}; b.lin_constraint_names_pretty = remove_underscores(b.lin_constraint_names); +% objectives +b.obj_names = {'LCOE','capex_design'}; +b.obj_names_pretty = {'LCOE','C_{design}'}; + +% indices [~,idxs_sort] = sort(b.var_names(1:end-1)); % alphabetical design variable indices idxs_recover = zeros(size(idxs_sort)); idxs_recover(idxs_sort) = 1:length(idxs_sort); % indices to recover unsorted variabes from sorted ones b.idxs_sort = idxs_sort; b.idxs_recover = idxs_recover; +% uuid b.filename_uuid = ''; % string to append to generated filenames to prevent parallel overlap +% calibrations of nominal values b.F_max_nom = find_nominal_inputs(b, mode, false); -b.X_noms(5) = b.F_max_nom; +b.X_noms(strcmp(b.var_names,'F_max')) = b.F_max_nom; -b.obj_names = {'LCOE','capex_design'}; -b.obj_names_pretty = {'LCOE','C_{design}'}; end \ No newline at end of file diff --git a/mdocean/optimization/find_nominal_inputs.m b/mdocean/optimization/find_nominal_inputs.m index e5280491..2c01528e 100644 --- a/mdocean/optimization/find_nominal_inputs.m +++ b/mdocean/optimization/find_nominal_inputs.m @@ -13,16 +13,17 @@ if display_on % check feasibility X = [b.X_noms; 1]; - X(5) = F_max_nom; + idx_F = strcmp(b.var_names,'F_max'); + X(idx_F) = F_max_nom; [LCOE, P_var, ~, g] = simulation(X,p) [feasible,failed] = is_feasible(g, X, p, b) % display x output - array2table(F_max_nom,'VariableNames',{'F_max (1e6 N)','B_p (1e6 Ns/m)','w_n (rad/s)'}) + array2table(F_max_nom,'VariableNames',{'F_max (1e6 N)'}) % display y output [~,y] = errFunc(x,y_desired,p,b); - results = round([y; y_desired] ./ [1000 1000 1],3,'significant'); + results = round([y; y_desired] ./ 1000,3,'significant'); array2table(results,'RowNames',{'Sim Output','RM3 Actual'},... 'VariableNames','Max Powertrain Force Error (-)') end @@ -31,7 +32,10 @@ function err = errFunc(F_max_in,p,b) X = [b.X_noms; 1]; - X(5) = F_max_in; + idx_F = strcmp(b.var_names,'F_max'); + X(idx_F) = F_max_in; + [~, ~, ~, g] = simulation(X, p); - err = g(12)^2; + idx_F_constr = strcmp(b.constraint_names,'irrelevant_max_force'); + err = g(idx_F_constr)^2; end diff --git a/mdocean/optimization/gradient_optim.m b/mdocean/optimization/gradient_optim.m index bb24bab9..5efc1c8b 100644 --- a/mdocean/optimization/gradient_optim.m +++ b/mdocean/optimization/gradient_optim.m @@ -1,4 +1,4 @@ -function [Xs_opt, objs_opt, flags, probs] = gradient_optim(x0_input,p,b,which_objs) +function [Xs_opt, objs_opt, flags, probs, vals] = gradient_optim(x0_input,p,b,which_objs) if nargin == 0 % set default parameters if function is run without input @@ -24,6 +24,8 @@ % create optimization variables for each of the design variables sz = [1 1]; % create scalar variables +% +assert(length(b.X_mins)==12) % this code currently assumes hardcoded number of design 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)); @@ -36,8 +38,6 @@ 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',... @@ -52,15 +52,15 @@ % iterate through material choices for matl = 1%1:2:3 %b.M_min : b.M_max - 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); + X = [x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 matl]; + [Xs_opt, objs_opt, flags, probs, vals] = optimize_both_objectives(X,p,b,x0_input,opts,ploton,which_objs); end end %% -function [Xs_opt, objs_opt, flags, probs] = optimize_both_objectives(X,p,b,x0_input,opts,ploton,which_objs) +function [Xs_opt, objs_opt, flags, probs, vals] = optimize_both_objectives(X,p,b,x0_input,opts,ploton,which_objs) num_constraints = length(b.constraint_names); num_objectives = length(which_objs); @@ -135,12 +135,17 @@ end 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 + [out(1),out(2),~,~,val] = simulation(X_opt,p); % rerun sim assert(out(which_obj) == obj_opt) % check correct reordering of X_opt elements Xs_opt(:,i) = X_opt; objs_opt(i) = obj_opt; flags(i) = flag; + if i==1 + vals = val; + else + vals(i) = val; + end % Post process if ploton diff --git a/mdocean/optimization/max_avg_power.m b/mdocean/optimization/max_avg_power.m new file mode 100644 index 00000000..7c7886be --- /dev/null +++ b/mdocean/optimization/max_avg_power.m @@ -0,0 +1,18 @@ +function [X_opt,val,flag] = max_avg_power(p,b) + +% hack that makes cost constant, so minimizing LCOE is actually maximizing average power +p_mod = p; +p_mod.cost_perN_mult = 0; +p_mod.cost_perW_mult = 0; +p_mod.cost_perkg_mult = [0 0 0]; + + +x0 = b.X_start_struct; + +% run LCOE minimization (effectively power maximization due to hack above) +[X_opt,~,flag] = gradient_optim(x0,p_mod,b,1); + +% plug back into unmodified simulation to get unsaturated power +[~, ~, ~, ~, val] = simulation(X_opt,p); + +end \ No newline at end of file diff --git a/mdocean/optimization/multiobjective/pareto_curve_heuristics.m b/mdocean/optimization/multiobjective/pareto_curve_heuristics.m index a82e5566..d3fc3c9a 100644 --- a/mdocean/optimization/multiobjective/pareto_curve_heuristics.m +++ b/mdocean/optimization/multiobjective/pareto_curve_heuristics.m @@ -1,5 +1,6 @@ function pareto_curve_heuristics() - p = parameters(); + close all + p0 = parameters(); b = var_bounds(); p_w = parameters('wecsim'); b_w = var_bounds('wecsim'); @@ -8,48 +9,68 @@ function pareto_curve_heuristics() d=dir("**/pareto_search_results*"); load(d(end).name) + if ~isequaln(p,p0) + warning(['You are loading results with different parameters than your ' ... + 'local machine right now. WecSim validation results (p_w) may be incorrect.']) + end + if ~exist('tol','var') tol = 1e-6; end - constraint_active_plot(residuals,fval,tol,b) + + new_objs = true; % switch between LCOE-Pvar and capex-Pavg + + constraint_active_plot(residuals,fval,tol,b,new_objs); 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 LCOE = fval(:,1); Pvar = fval(:,2); - - new_objs = true; % switch between LCOE-Pvar and capex-Pavg - [J1, minJ1, idx_best_J1, J1_nom, J1_nom_sim, J1_solar, J1_balanced,... - J2, minJ2, idx_best_J2, J2_nom, J2_nom_sim, J2_solar, J2_balanced,... +% if new_objs +% [X_max_pwr,val] = max_avg_power(p,b); +% [LCOE_max_pwr,P_var_max_pwr] = simulation(X_max_pwr,p); % fixme do I need a constraint to ensure P_max is tight, like I have for F_max? +% LCOE(end+1) = LCOE_max_pwr; +% Pvar(end+1) = P_var_max_pwr; +% X(end+1,:) = X_max_pwr; +% end + + [J1, bestJ1, idx_best_J1, J1_nom, J1_nom_sim, J1_solar, J1_balanced,... + J2, bestJ2, idx_best_J2, J2_nom, J2_nom_sim, J2_solar, J2_balanced,... x_best_J1, x_best_J2, x_nom, x_balanced, idxo] = process_pareto_front(LCOE,Pvar,X,p,p_w,b,b_w, new_objs); %% super simple "pareto" plot of just single objective optimizations showSingleObj = true; showImages = false; - pareto_plot(J1, minJ1, idx_best_J1, J1_nom, J1_nom_sim, J1_solar, NaN*J1_balanced,... - J2, minJ2, idx_best_J2, J2_nom, J2_nom_sim, J2_solar, NaN*J2_balanced,... - x_best_J1, x_best_J2, x_nom, x_balanced, [], showSingleObj, showImages, p, new_objs) + showLCOEContours = false; + pareto_plot(J1, bestJ1, idx_best_J1, J1_nom, J1_nom_sim, J1_solar, NaN*J1_balanced,... + J2, bestJ2, idx_best_J2, J2_nom, J2_nom_sim, J2_solar, NaN*J2_balanced,... + x_best_J1, x_best_J2, x_nom, x_balanced, [], showSingleObj, ... + showImages, showLCOEContours, p, new_objs) %% simple pareto plot showSingleObj = false; showImages = false; - pareto_plot(J1, minJ1, idx_best_J1, J1_nom, J1_nom_sim, J1_solar, J1_balanced,... - J2, minJ2, idx_best_J2, J2_nom, J2_nom_sim, J2_solar, J2_balanced,... - x_best_J1, x_best_J2, x_nom, x_balanced, idxo, showSingleObj, showImages, p, new_objs) + showLCOEContours = true; + pareto_plot(J1, bestJ1, idx_best_J1, J1_nom, J1_nom_sim.*[1 NaN], J1_solar, J1_balanced,... + J2, bestJ2, idx_best_J2, J2_nom, J2_nom_sim.*[1 NaN], J2_solar, J2_balanced,... + x_best_J1, x_best_J2, x_nom, x_balanced, idxo, showSingleObj, ... + showImages, showLCOEContours, p, new_objs) %% plot pareto front with annotations and embedded images of three recommended designs showSingleObj = true; showImages = true; - pareto_plot(J1, minJ1, idx_best_J1, J1_nom, J1_nom_sim, J1_solar, J1_balanced,... - J2, minJ2, idx_best_J2, J2_nom, J2_nom_sim, J2_solar, J2_balanced,... - x_best_J1, x_best_J2, x_nom, x_balanced, idxo, showSingleObj, showImages, p, new_objs) + showLCOEContours = false; + pareto_plot(J1, bestJ1, idx_best_J1, J1_nom, J1_nom_sim.*[1 NaN], J1_solar, J1_balanced,... + J2, bestJ2, idx_best_J2, J2_nom, J2_nom_sim.*[1 NaN], J2_solar, J2_balanced,... + x_best_J1, x_best_J2, x_nom, x_balanced, idxo, showSingleObj, ... + showImages, showLCOEContours, p, new_objs) %% plots for DVs as a fn of percent along the pareto J1_max = Inf;%p.LCOE_max; - design_heuristics_plot(J1, minJ1, idx_best_J1, x_best_J1, ... - J2, minJ2, idx_best_J2, X, idxo, J1_max, b.var_names_pretty(1:end-1),new_objs) + design_heuristics_plot(J1, bestJ1, idx_best_J1, x_best_J1, ... + J2, bestJ2, idx_best_J2, X, idxo, J1_max, b.var_names_pretty(1:end-1),new_objs) end %% @@ -74,7 +95,7 @@ function pareto_curve_heuristics() capex_design = Pvar; % second sim output means something different J2 = capex_design; - J2_fieldname = 'capex_design'; + J2_fieldname = 'J_capex_design'; J2_fieldidx = 4; J2_mult = 1; % $M @@ -169,36 +190,39 @@ function pareto_curve_heuristics() end %% -function [] = pareto_plot(LCOE,minLCOE,idx_best_LCOE,LCOE_nom, LCOE_nom_sim, LCOE_solar, LCOE_balanced,... - Pvar,minPvar,idx_best_Pvar,P_var_nom,P_var_nom_sim,P_var_solar,P_var_balanced,... - x_best_LCOE,x_best_Pvar,x_nom,x_balanced,idxo,showSingleObj,showImages,p,new_objs) +function [] = pareto_plot(J1, bestJ1, idx_best_J1, J1_nom, J1_nom_sim, J1_solar, J1_balanced,... + J2, bestJ2, idx_best_J2, J2_nom, J2_nom_sim, J2_solar, J2_balanced,... + x_best_J1, x_best_J2, x_nom, x_balanced, idxo, showSingleObj,... + showImages, showLCOEContours, p, new_objs) figure % overall pareto front - plot(LCOE(idxo),Pvar(idxo),'bs','MarkerFaceColor','b','HandleVisibility','off') + plot(J1(idxo),J2(idxo),'bs','MarkerFaceColor','b','HandleVisibility','off') hold on % utopia point - plot(minLCOE,minPvar,'gp','MarkerFaceColor','g','MarkerSize',20,'HandleVisibility','off') + plot(bestJ1,bestJ2,'gp','MarkerFaceColor','g','MarkerSize',20,'HandleVisibility','off') % RM3 nominal reference - report - plot(LCOE_nom(1),P_var_nom(1),'rd','HandleVisibility','off') - plot(LCOE_nom_sim(1),P_var_nom_sim(1),'rs','HandleVisibility','off') + plot(J1_nom(1), J2_nom(1), 'rd','HandleVisibility','off') + plot(J1_nom_sim(1),J2_nom_sim(1),'rs','HandleVisibility','off') % RM3 nominal reference - wecsim - plot(LCOE_nom(2),P_var_nom(2),'md','HandleVisibility','off') - plot(LCOE_nom_sim(2),P_var_nom_sim(2),'ms','HandleVisibility','off') + plot(J1_nom(2), J2_nom(2), 'md','HandleVisibility','off') + plot(J1_nom_sim(2),J2_nom_sim(2),'ms','HandleVisibility','off') if showSingleObj % black squares for 3 ref points - 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') + plot(bestJ1,J2(idx_best_J1),'ks','HandleVisibility','off') + plot(J1(idx_best_J2),bestJ2,'ks','HandleVisibility','off') + plot(J1_balanced,J2_balanced,'ks','HandleVisibility','off') end % axis labels if new_objs xlabel('Average Electrical Power (kW)') ylabel('Structural and PTO Cost ($M)') + xlim([80 300]) + ylim([.5 3]) else xlabel('LCOE ($/kWh)') ylabel('Power Variation (%)') @@ -211,105 +235,104 @@ function pareto_curve_heuristics() % solar reference yellow = [1, .87, .2]; - plot(LCOE_solar, P_var_solar,'o','MarkerSize',12,'MarkerEdgeColor',yellow,... + plot(J1_solar, J2_solar,'o','MarkerSize',12,'MarkerEdgeColor',yellow,... 'MarkerFaceColor',yellow,'HandleVisibility','off') % for the yellow color to work, do not use improvePlot below here % text labels sz = 14; - text(LCOE_nom(1)+.03,P_var_nom(1),'RM3 Report','FontSize',sz) - text(LCOE_nom(1)+.01,P_var_nom(1)-5,'Actual [10]','FontSize',sz) - text(LCOE_nom_sim(1)-.02,P_var_nom_sim(1)+5,'RM3 Report','FontSize',sz) - text(LCOE_nom_sim(1)+.03,P_var_nom_sim(1),'Sim','FontSize',sz) - - text(LCOE_nom(2)+.03,P_var_nom(2),'RM3 WecSim','FontSize',sz) - text(LCOE_nom(2)+.01,P_var_nom(2)-5,'Actual','FontSize',sz) - text(LCOE_nom_sim(2)-.02,P_var_nom_sim(2)+5,'RM3 WecSim','FontSize',sz) - text(LCOE_nom_sim(2)+.03,P_var_nom_sim(2),'Sim','FontSize',sz) - - text(minLCOE+.03,minPvar-2,'Utopia Point','FontSize',sz) - text(LCOE_solar+.03,P_var_solar,'Solar','FontSize',sz) + text(J1_nom(1)+10,J2_nom(1),'RM3 Report','FontSize',sz) + text(J1_nom(1)+.01,J2_nom(1)-5,'Actual [10]','FontSize',sz) + text(J1_nom_sim(1)-.02,J2_nom_sim(1)+5,'RM3 Report','FontSize',sz) + text(J1_nom_sim(1)+.03,J2_nom_sim(1),'Sim','FontSize',sz) + + text(J1_nom(2)+.03,J2_nom(2),'RM3 WecSim','FontSize',sz) + text(J1_nom(2)+.01,J2_nom(2)-5,'Actual','FontSize',sz) + text(J1_nom_sim(2)-.02,J2_nom_sim(2)+5,'RM3 WecSim','FontSize',sz) + text(J1_nom_sim(2)+.03,J2_nom_sim(2),'Sim','FontSize',sz) + + text(bestJ1+.03,bestJ2-2,'Utopia Point','FontSize',sz) + text(J1_solar+.03,J2_solar,'Solar','FontSize',sz) if showSingleObj - text(LCOE(idx_best_LCOE)+.03,Pvar(idx_best_LCOE),'Cheapest','FontSize',sz) - 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) + text(J1(idx_best_J1)+.03,J2(idx_best_J1),'Cheapest','FontSize',sz) + text(J1(idx_best_J2)+.03,J2(idx_best_J2)-3,'Least Variable','FontSize',sz) + text(J1_balanced-.15,J2_balanced+5,'Balanced Design','FontSize',sz) end - if showImages - mini_plot_size = [.2 .22]; - % small corner pictures of best geometries - % upper left - axes('Position',[.28 .6 mini_plot_size]) - box on - visualize_geometry(x_best_LCOE,p,true); - set(gca,'XTickLabel',[],'YTickLabel',[]) - - % lower right - axes('Position',[.51 .23 mini_plot_size]) - box on - visualize_geometry(x_best_Pvar,p,true); - set(gca,'XTickLabel',[],'YTickLabel',[]) + if showLCOEContours + LCOE_min = 0.272; + LCOE_nom = 0.76; % fixme hardcoded + overlay_LCOE(p, LCOE_nom, LCOE_min) + end + + if showImages % small pictures of best geometries + + upper_left = [.28 .6]; + mini_plot(upper_left,x_best_J1,p) + + lower_right = [.51 .23]; + mini_plot(lower_right,x_best_J2,p) - % balanced - axes('Position',[.10 .28 mini_plot_size]) - box on - visualize_geometry(x_balanced,p,true); - set(gca,'XTickLabel',[],'YTickLabel',[]) + balanced_pos = [.10 .28]; + mini_plot(balanced_pos,x_balanced,p) + + report_pos = [.7 .53]; + mini_plot(report_pos,x_nom(:,1),p) - % RM3 report - axes('Position',[.7 .53 mini_plot_size]) - box on - visualize_geometry(x_nom(:,1),p,true); - set(gca,'XTickLabel',[],'YTickLabel',[]) - - % RM3 WecSim - axes('Position',[.8 .7 mini_plot_size]) - box on - visualize_geometry(x_nom(:,2),p,true); - set(gca,'XTickLabel',[],'YTickLabel',[]) + wecsim_pos = [.8 .7]; + mini_plot(wecsim_pos,x_nom(:,2),p) end + +end + +function mini_plot(pos, x, p) + mini_plot_size = [.2 .22]; + axes('Position',[pos mini_plot_size]) + box on + visualize_geometry(x,p,true); + set(gca,'XTickLabel',[],'YTickLabel',[]) end %% -function [] = design_heuristics_plot(overallLCOE, minLCOE, idx_best_LCOE, x_best_LCOE, ... - overallPvar, minPvar, idx_best_Pvar, ... - overallX, idxo, LCOE_max, var_names, new_objs) - - LCOE_pareto = overallLCOE(idxo); - Pvar_pareto = overallPvar(idxo); - [LCOE_pareto_sorted,idx_sort] = sort(LCOE_pareto(LCOE_pareto0),w_below_pk(P_below_pk>0),P_half); + w_half_above = interp1(P_above_pk(P_above_pk>0),w_above_pk(P_above_pk>0),P_half); + delta_w = w_half_above - w_half_below; + +% figure +% plot(w,P) +% hold on +% plot(w(idx_max),max_power,'mp') +% plot(w_half_below,P_half,'ro') +% plot(w_half_above,P_half,'go') + +end + +function b = fix_constraints(p,b) + desired_length = 20 + numel(p.JPD) + length(p.T_struct); + len = length(b.constraint_names); + if len < desired_length + b.constraint_names(end+1 : end+desired_length-len) = 'slamming'; + b.constraint_names_pretty = remove_underscores(b.constraint_names); + elseif len > desired_length + b.constraint_names = b.constraint_names(1:desired_length); + b.constraint_names_pretty = remove_underscores(b.constraint_names); + end end \ No newline at end of file diff --git a/mdocean/plots/all_figures.m b/mdocean/plots/all_figures.m index f8e56ca7..54e06a34 100644 --- a/mdocean/plots/all_figures.m +++ b/mdocean/plots/all_figures.m @@ -7,8 +7,12 @@ % since prob2struct needs unique filenames for code generation end +date = datestr(now,'yyyy-mm-dd_HH.MM.SS'); +save_folder = ['../test-results/' date '/']; +mkdir(save_folder) + num_figs = 29; -num_tabs = 7; +num_tabs = 8; fig_names = cell([1,num_figs]); tab_names = cell([1,num_tabs]); success_criterion = {}; @@ -19,9 +23,14 @@ which_tabs = 1:num_tabs; end -fig_output = gobjects(1, length(which_figs)); +if isempty(which_figs) + which_figs = 0; +end +if isempty(which_tabs) + which_tabs = 0; +end -%fig_output = cell(1, length(which_figs)); +fig_output = gobjects(1, length(which_figs)); tab_output(1, 1:length(which_tabs)) = {table()}; @@ -90,27 +99,19 @@ fig_names{6} = 'Fig. 6: hydro coeffs vs freq'; if any(which_figs == 6) fig6 = figure; - % fixme - not implemented + hydro_coeff_err() fig_output(which_figs==6) = fig6; end -%% figure 7 - drag DF +%% figure 7, 8 - drag DF, saturation time signal fig_names{7} = 'Fig. 7: drag describing function'; -if any(which_figs == 7) - fig7 = figure; - % fixme - not implemented - fig_output(which_figs==7) = fig7; -end - -%% figure 8, XX - saturation time signal, saturation alpha fig_names{8} = 'Fig. 8: force saturation time signal'; -%fig_names{4} = 'Fig. XX: saturation alpha'; -if any(which_figs == 8) % | which_figs == 4) +if any(which_figs == 7 | which_figs == 8) sin_desc_fcn_demo() - figX = gcf; - fig8 = figure(figX.Number-1); + fig7 = gcf; + fig8 = figure(fig7.Number-2); + fig_output(which_figs==7) = fig7; fig_output(which_figs==8) = fig8; - %fig_output(which_figs==4) = figX; end %% figure 9 - JPD multiplication @@ -119,7 +120,7 @@ p = parameters(); b = var_bounds(); X = [b.X_noms; 1]; - plot_power_matrix(X,p) + plot_power_matrix(X,p,filename_uuid) fig9 = gcf; fig_output(which_figs==9) = fig9; end @@ -127,8 +128,8 @@ %% figure 10 - WecSim error breakdown fig_names{10} = 'Fig. 10: WECSim error breakdown'; if any(which_figs == 10) - % fixme - not implemented fig10 = figure; + imshow(imread("error_breakdown.png"),'Parent',axes(fig10)); fig_output(which_figs==10) = fig10; end @@ -141,6 +142,7 @@ fig_output(which_figs==11) = fig11; end + %% figure 12 - cost vs N WEC fig_names{12} = 'Fig. 12: cost vs N WEC'; if any(which_figs == 12) @@ -152,8 +154,8 @@ %% figure 13 - sim runtime fig_names{13} = 'Fig. 13: sim runtime'; if any(which_figs == 13) - % fixme - not implemented fig13 = figure; + imshow(imread("tree_map.png"),'Parent',axes(fig13)); fig_output(which_figs==13) = fig13; end @@ -189,23 +191,24 @@ end %% figure 22, 23, 24, 25, 26 - pareto front, design heuristics -fig_names{22} = 'Fig. 22: pareto front'; +fig_names{22} = 'Fig. 22: pareto front with design images'; fig_names{23} = 'Fig. 23: design heuristics'; fig_names{24} = 'Fig. 24: objective heuristics'; fig_names{25} = 'Fig. 25: pareto front with LCOE contours'; fig_names{26} = 'Fig. 26: constraint activity'; -if any(which_figs == 22 | which_figs == 23 | which_figs == 24) +if any(which_figs == 22 | which_figs == 23 | which_figs == 24 | which_figs == 25 | which_figs == 26) pareto_search(filename_uuid); pareto_curve_heuristics() fig24 = gcf; fig23 = figure(fig24.Number - 1); fig22 = figure(fig24.Number - 3); + fig25 = figure(fig24.Number - 4); fig26 = figure(fig24.Number - 6); % constraint activity fig26.Position = [1 41 1536 844.8000]; fig_output(which_figs==22) = fig22; fig_output(which_figs==23) = fig23; fig_output(which_figs==24) = fig24; - % fixme: 25 not implemented + fig_output(which_figs==25) = fig25; fig_output(which_figs==26) = fig26; end @@ -213,8 +216,8 @@ fig_names{27} = 'Fig. 27: overlaid geometry'; fig_names{28} = 'Fig. 28: overlaid hydro coeffs'; fig_names{29} = 'Fig. 29: probability CDF'; -if any(which_figs == 27 | which_figs == 28 | which_figs == 29 | which_tabs == 5) - tab5 = compare(filename_uuid); +if any(which_figs == 27 | which_figs == 28 | which_figs == 29 | which_tabs == 5 | which_tabs == 6) + [tab5,tab6] = compare(filename_uuid); n = gcf().Number; fig29 = figure(n-1); fig27 = figure(n-2); @@ -226,18 +229,29 @@ % fixme: 30 to 46 not implemented -%% table 1 - design variables table -tab_names{1} = 'Tab. 1: design variables'; +%% table 12 - validation table +tab_names{1} = 'Tab. 12: validation against nominal'; if any(which_tabs == 1) - b = var_bounds(); - tab1 = array2table([b.X_mins b.X_noms b.X_maxs], ... - 'VariableNames',{'Mins','Noms','Maxs'}, 'RowNames', b.var_names(1:end-1)); + [~,~,~,~,tab1a] = validate_nominal_RM3('report'); + display(tab1a) + [~,~,~,~,tab1b] = validate_nominal_RM3('wecsim'); + display(tab1b) + + % merge table 1a and 1b while preserving row order + sharedCols = intersect(tab1a.Properties.VariableNames, tab1b.Properties.VariableNames); + tab1a.RowNum = (1:length(tab1a.Properties.RowNames))'; + tab1b.RowNum = length(tab1a.Properties.RowNames) + (1:length(tab1b.Properties.RowNames))'; + + tab1 = outerjoin(tab1a, tab1b, 'Keys', [{'RowNum'},sharedCols], 'MergeKeys', true); + tab1 = removevars(tab1,'RowNum'); + tab1.Properties.RowNames = [tab1a.Properties.RowNames; tab1b.Properties.RowNames]; display(tab1) + tab_output{which_tabs==1} = tab1; end -%% table 2 - constraints table -tab_names{2} = 'Tab. 2: constraints'; +%% table 15 - constraints table +tab_names{2} = 'Tab. 15: constraints'; if any(which_tabs == 2) b = var_bounds(); tab2 = b.constraint_names'; @@ -245,58 +259,65 @@ tab_output{which_tabs==2} = tab2; end -%% table 3 - parameters table -tab_names{3} = 'Tab. 3: parameters'; +%% table 16 - design variables table +tab_names{3} = 'Tab. 16: design variables'; if any(which_tabs == 3) - [~,tab3] = parameters(); + b = var_bounds(); + tab3 = array2table([b.X_mins b.X_noms b.X_maxs], ... + 'VariableNames',{'Mins','Noms','Maxs'}, 'RowNames', b.var_names(1:end-1)); display(tab3) tab_output{which_tabs==3} = tab3; end -%% table 4 - validation table -tab_names{4} = 'Tab. 4: validation against nominal'; +%% table 17 - parameters table +tab_names{4} = 'Tab. 17: parameters'; if any(which_tabs == 4) - [~,~,~,~,tab4a] = validate_nominal_RM3('report'); - display(tab4a) - [~,~,~,~,tab4b] = validate_nominal_RM3('wecsim'); - display(tab4b) - - % merge table 4a and 4b while preserving row order - sharedCols = intersect(tab4a.Properties.VariableNames, tab4b.Properties.VariableNames); - tab4a.RowNum = (1:length(tab4a.Properties.RowNames))'; - tab4b.RowNum = length(tab4a.Properties.RowNames) + (1:length(tab4b.Properties.RowNames))'; - - tab4 = outerjoin(tab4a, tab4b, 'Keys', [{'RowNum'},sharedCols], 'MergeKeys', true); - tab4 = removevars(tab4,'RowNum'); - tab4.Properties.RowNames = [tab4a.Properties.RowNames; tab4b.Properties.RowNames]; + [~,tab4] = parameters(); display(tab4) - tab_output{which_tabs==4} = tab4; end -%% table 5 - optimal DVs for 4 designs -tab_names{5} = 'Tab. 5: optimal DVs for 4 designs'; +%% table 19 - optimal DVs for various designs +tab_names{5} = 'Tab. 19: optimal DVs for various designs'; if any(which_tabs == 5) - % computation above with figures 9-10 + % computation above with figures 27-29 display(tab5); tab_output{which_tabs==5} = tab5; + table2latex(tab5,[save_folder 'table_19.tex']) end -%% table 6 - optimal DVs for 4 locations -tab_names{6} = 'Tab. 6: optimal DVs for 4 locations'; +%% table 20 - optimal outputs for various designs +tab_names{6} = 'Tab. 20: optimal outputs for various designs'; if any(which_tabs == 6) - tab6 = location_sensitivity(filename_uuid); + % computation above with figures 27-29 display(tab6); tab_output{which_tabs==6} = tab6; - location_flags = tab6(strcmp(tab6.Row,'flag'),:).Variables; - success_criterion(end+1) = {location_flags}; + table2latex(tab6,[save_folder 'table_20.tex']) end -%% paragraph 4.2 - convergence for different x0s -tab_names{7} = 'Tab. 7: convergence for different x0s'; +%% table 21 - convergence for different x0s +tab_names{7} = 'Tab. 21: convergence for different x0s'; if any(which_tabs == 7) tab7 = gradient_mult_x0(filename_uuid); tab_output{which_tabs==7} = tab7; + table2latex(tab7,[save_folder 'table_21.tex']) end +%% table 22 - optimal DVs for 4 locations +tab_names{8} = 'Tab. 22: optimal DVs for 4 locations'; +if any(which_tabs == 8) + tab8 = location_sensitivity(filename_uuid); + display(tab8); + tab_output{which_tabs==8} = tab8; + location_flags = str2double(tab8(strcmp(tab8.Row,'flag'),:).Variables); + success_criterion(end+1) = {location_flags}; + + idx_remove = ismember(tab8.Row,{'flag','Optimal Material index'}); + tab8latex = tab8(~idx_remove,:); + colspec = '>{\centering\arraybackslash}p{0.20\linewidth}>{\centering\arraybackslash}p{0.08\linewidth}>{\centering\arraybackslash}p{0.15\linewidth}>{\centering\arraybackslash}p{0.15\linewidth}>{\centering\arraybackslash}p{0.15\linewidth}>{\centering\arraybackslash}p{0.18\linewidth}'; + firstrow = '&& \multicolumn{4}{c}{Location}\\ \cline{3-6}'; + table2latex(tab8latex,[save_folder 'table_22.tex'],colspec,firstrow) +end + + end \ No newline at end of file diff --git a/mdocean/plots/compare.m b/mdocean/plots/compare.m index ffadc3d0..27742dce 100644 --- a/mdocean/plots/compare.m +++ b/mdocean/plots/compare.m @@ -1,67 +1,85 @@ -function [DV_table] = compare(filename_uuid) +function [DV_table,out_table] = compare(filename_uuid) + +if nargin==0 + filename_uuid = ''; +end p = parameters(); b = var_bounds(); b.filename_uuid = filename_uuid; x0_input = b.X_start_struct; -Xs_opt = gradient_optim(x0_input,p,b); -cheapest = Xs_opt(:,1); -minVariation = Xs_opt(:,2); +[Xs_opt,~,~,~,vals_opt] = gradient_optim(x0_input,p,b); +X_minLCOE = Xs_opt(:,1); +X_minCapex = Xs_opt(:,2); +val_minLCOE = vals_opt(1); +val_minCapex = vals_opt(2); + +[X_maxPower,val_maxPower] = max_avg_power(p,b); p.LCOE_max = 0.1; -balanced = gradient_optim(x0_input,p,b,2); +[X_balanced,~,~,~,val_balanced] = gradient_optim(x0_input,p,b,2); + +X_nom = [b.X_noms' 1]; +[~,~,~,~,val_nom] = simulation(X_nom,p); -nominal = [b.X_noms' 1]; %% -X = [nominal;cheapest';minVariation';balanced']; +X = [X_nom; X_minLCOE'; X_minCapex'; X_maxPower'; X_balanced']; +vals = [val_nom, val_minLCOE, val_minCapex, val_maxPower, val_balanced]; -col = {'k','b','r','g'}; +titles = {'Nominal','Min LCOE','Min Capex','Max Power','Balanced'}; +color = {'k','b','r','g','m'}; +num_designs = length(titles); %% geometry comparison figure -for i=1:4 +for i=1:num_designs x = X(i,:); hold on - visualize_geometry(x,p,false,col{i}) + visualize_geometry(x,p,false,color{i}) end -legend('Nominal','Min LCOE','Min Variation','Balanced') +legend(titles) %% power probability comparison figure t = tiledlayout(2,1); t.TileSpacing = 'compact'; ylabel(t,'Probability (-)','FontWeight','bold') -for i=1:4 +for i=1:num_designs x = X(i,:); hold on - t = power_PDF(x,p,col{i},t); + t = power_PDF(x,p,color{i},t); end -legend('Nominal','Min LCOE','Min Variation','Balanced','location','best') +legend(titles{:},'location','best') %% power matrix comparison -titles = {'Nominal','Min LCOE','Min Variation','Balanced'}; + figure -for i=1:4 +for i=1:num_designs x = X(i,:); [~,~,P_matrix] = simulation(x, p); P_matrix = P_matrix / 1e3; % convert W to kW [T,Hs] = meshgrid(p.T,p.Hs); - subplot(2,2,i); + subplot(2,3,i); hold on - contourf(T,Hs,P_matrix .* p.JPD/100,[0 5 10:10:60]) + contourf(T,Hs,P_matrix .* p.JPD/100) improvePlot grid on title(titles{i}) - caxis([0 60]) + caxis([0 15]) ax = gca; - if i==2 || i==4 + second_col = i==2 || i==5; + third_col = i==3; + if second_col xoffset = -ax.Position(3)*.25; + elseif third_col + xoffset = -ax.Position(3)*.5; else xoffset = 0; end - if i==3 || i==4 + second_row = i>3; + if second_row yoffset = ax.Position(4)*.2; else yoffset = 0; @@ -79,6 +97,13 @@ %% design variable table DV_table = array2table(X.', ... - 'VariableNames',titles, 'RowNames', b.var_names); + 'VariableNames',titles, 'RowNames', b.var_names_pretty); + +%% output table +temp_table = struct2table(vals,'RowNames',titles); +scalar_vals = varfun(@isnumeric, temp_table, 'OutputFormat', 'uniform'); +out_table = rows2vars(temp_table(:,scalar_vals)); +out_table.Properties.RowNames = out_table.OriginalVariableNames; +out_table = removevars(out_table,'OriginalVariableNames'); end \ No newline at end of file diff --git a/mdocean/plots/constraint_active_plot.m b/mdocean/plots/constraint_active_plot.m index 0729fed4..5012ba27 100644 --- a/mdocean/plots/constraint_active_plot.m +++ b/mdocean/plots/constraint_active_plot.m @@ -1,4 +1,9 @@ -function [idx] = constraint_active_plot(residuals,fval,tol,b) +function [idx] = constraint_active_plot(residuals,fval,tol,b,reversed) + + if nargin<5 + reversed = false; + end + lb_active = abs(residuals.lower) < tol; ub_active = abs(residuals.upper) < tol; nlcon_active = abs(residuals.ineqnonlin) < tol; @@ -10,8 +15,15 @@ idx_slamming_after = idx_slamming & ~idx_slamming_first; nlcon_active(:,idx_slamming_first) = any(nlcon_active(:,idx_slamming),2); nlcon_active(:,idx_slamming_after) = []; + constraint_names_mod = b.constraint_names_pretty(~idx_slamming_after); + constraint_names_mod(idx_slamming_first) = {'Prevent Slamming'}; - [~,idx] = sort(fval(:,1)); % order by increasing LCOE + if reversed + dir = 'descend'; + else + dir = 'ascend'; + end + [~,idx] = sort(fval(:,1),dir); % order by increasing LCOE (decreasing if reversed) figure tiledlayout(2,2); @@ -34,7 +46,7 @@ spy(nlcon_active(idx,:)') title('Nonlinear Constraint Active') grid on - set(gca,'ytick',1:size(nlcon_active,2),'yticklabel',b.constraint_names_pretty(~idx_slamming_after)) + set(gca,'ytick',1:size(nlcon_active,2),'yticklabel',constraint_names_mod) set(gca, 'PlotBoxAspectRatio', [1.5 1 1]) nexttile diff --git a/mdocean/plots/non_matlab_figs/dimensions.pdf b/mdocean/plots/non_matlab_figs/dimensions.pdf index f2c7000a..2ef7ac5a 100644 Binary files a/mdocean/plots/non_matlab_figs/dimensions.pdf and b/mdocean/plots/non_matlab_figs/dimensions.pdf differ diff --git a/mdocean/plots/non_matlab_figs/error_breakdown.png b/mdocean/plots/non_matlab_figs/error_breakdown.png new file mode 100644 index 00000000..1dbc1d93 Binary files /dev/null and b/mdocean/plots/non_matlab_figs/error_breakdown.png differ diff --git a/mdocean/plots/non_matlab_figs/tree_map.png b/mdocean/plots/non_matlab_figs/tree_map.png new file mode 100644 index 00000000..ba85e8df Binary files /dev/null and b/mdocean/plots/non_matlab_figs/tree_map.png differ diff --git a/mdocean/plots/non_matlab_figs/xdsm.pdf b/mdocean/plots/non_matlab_figs/xdsm.pdf index ca542545..b6f4dda1 100644 Binary files a/mdocean/plots/non_matlab_figs/xdsm.pdf and b/mdocean/plots/non_matlab_figs/xdsm.pdf differ diff --git a/mdocean/plots/non_matlab_figs/xdsm.tikz b/mdocean/plots/non_matlab_figs/xdsm.tikz index 2acc7669..c3eeecf8 100644 --- a/mdocean/plots/non_matlab_figs/xdsm.tikz +++ b/mdocean/plots/non_matlab_figs/xdsm.tikz @@ -22,6 +22,7 @@ \node [Optimization] (opt) {$\text{Optimizer}$};& \node [DataInter] (opt-geom) {$\begin{array}{c}\text{Dimensions,} \\ \text{Thicknesses}\end{array}$};& \node [DataInter] (opt-hydro) {$\text{Dimensions}$};& +& \node [DataInter] (opt-dynam) {$\begin{array}{c}\text{Generator} \\ \text{ratings}\end{array}$};& \node [DataInter] (opt-struct) {$\begin{array}{c}\text{Dimensions,} \\ \text{thicknesses}\end{array}$};& \node [DataInter] (opt-econ) {$\begin{array}{c}\text{Generator} \\ \text{ratings}\end{array}$};& @@ -32,6 +33,7 @@ & \node [Function] (geom) {$\text{Geometry}$};& & +& \node [DataInter] (geom-dynam) {$\text{Mass}$};& & \node [DataInter] (geom-econ) {$\text{Material volume}$};& @@ -42,6 +44,7 @@ & & \node [Function] (hydro) {$\text{Hydrodynamics}$};& +& \node [DataInter] (hydro-dynam) {$\begin{array}{c}\text{Hydrodynamic} \\ \text{coefficients}\end{array}$};& & & @@ -52,12 +55,25 @@ & & & +\node [MDA] (solver) {$\text{Iteration}$};& +\node [DataInter] (solver-dynam) {$\begin{array}{c}\text{Dynamic} \\ \text{response guess}\end{array}$};& +& +& +& +\\ +%Row 4 +& +& +& +& +\node [DataInter] (dynam-solver) {$\begin{array}{c}\text{Dynamic} \\ \text{response residual}\end{array}$};& \node [Function] (dynam) {$\begin{array}{c}\text{Dynamics} \\ \text{and Control}\end{array}$};& \node [DataInter] (dynam-struct) {$\text{Loads}$};& \node [DataInter] (dynam-econ) {$\text{Power}$};& \node [DataInter] (dynam-F) {$\text{Power}$};& \node [DataInter] (dynam-G) {$\text{Amplitude constraints}$};\\ -%Row 4 +%Row 5 +& & & & @@ -67,7 +83,8 @@ & & \node [DataInter] (struct-G) {$\text{Structural constraints}$};\\ -%Row 5 +%Row 6 +& & & & @@ -77,7 +94,7 @@ \node [Function] (econ) {$\text{Economics}$};& \node [DataInter] (econ-F) {$\begin{array}{c}\text{LCOE,} \\ \text{Capital cost}\end{array}$};& \\ -%Row 6 +%Row 7 \node [DataIO] (left_output_F) {$J^*$};& \node [DataInter] (F-opt) {$J$};& & @@ -85,9 +102,10 @@ & & & -\node [Function] (F) {$Objective$};& +& +\node [Function] (F) {$\text{Objective}$};& \\ -%Row 7 +%Row 8 \node [DataIO] (left_output_G) {$g^*$};& \node [DataInter] (G-opt) {$g$};& & @@ -96,8 +114,10 @@ & & & -\node [Function] (G) {$Constraints$};\\ -%Row 8 +& +\node [Function] (G) {$\text{Constraints}$};\\ +%Row 9 +& & & & @@ -120,6 +140,8 @@ (opt) edge [DataLine] (opt-dynam) (opt) edge [DataLine] (opt-struct) (opt) edge [DataLine] (opt-econ) +(solver) edge [DataLine] (solver-dynam) +(dynam) edge [DataLine] (dynam-solver) (geom) edge [DataLine] (geom-dynam) (geom) edge [DataLine] (geom-econ) (hydro) edge [DataLine] (hydro-dynam) @@ -141,6 +163,8 @@ (opt-dynam) edge [DataLine] (dynam) (opt-struct) edge [DataLine] (struct) (opt-econ) edge [DataLine] (econ) +(solver-dynam) edge [DataLine] (dynam) +(dynam-solver) edge [DataLine] (solver) (geom-dynam) edge [DataLine] (dynam) (geom-econ) edge [DataLine] (econ) (hydro-dynam) edge [DataLine] (dynam) diff --git a/mdocean/plots/plot_power_matrix.m b/mdocean/plots/plot_power_matrix.m index 5ed129c3..879501fd 100644 --- a/mdocean/plots/plot_power_matrix.m +++ b/mdocean/plots/plot_power_matrix.m @@ -1,33 +1,69 @@ -function plot_power_matrix(X,p) +function plot_power_matrix(X,p,filename_uuid) [~,~,P_matrix] = simulation(X, p); P_matrix = P_matrix / 1e3; % convert W to kW +[CW_over_CW_max, P_wave, CW_max] = check_max_CW(filename_uuid); [T,Hs] = meshgrid(p.T,p.Hs); figure -subplot(1,3,1) -contourf(T,Hs,P_matrix); +% subplot(1,3,1) +% contourf(T,Hs,P_matrix); +% xlabel('Wave Period T (s)') +% ylabel('Wave Height Hs (m)') +% title('Power (kW)') +% colorbar +% +% subplot(1,3,2) +% contourf(T,Hs,p.JPD) +% xlabel('Wave Period T (s)') +% ylabel('Wave Height Hs (m)') +% title('Probability (%)') +% colorbar +% +% subplot(1,3,3) +% contourf(T,Hs,P_matrix .* p.JPD/100) +% xlabel('Wave Period T (s)') +% ylabel('Wave Height Hs (m)') +% title('Weighted Power (kW)') +% colorbar + +subplot(1,5,1) +contourf(T,Hs,P_wave); +xlabel('Wave Period T (s)') +ylabel('Wave Height Hs (m)') +title('Raw Wave Power Density (W/m)','FontSize',12) +colorbar + +subplot(1,5,2) +contourf(T,Hs,CW_over_CW_max) +xlabel('Wave Period T (s)') +ylabel('Wave Height Hs (m)') +title('Device Capture Efficiency','FontSize',12) +colorbar + +subplot(1,5,3) +contourf(T,Hs,CW_max) xlabel('Wave Period T (s)') ylabel('Wave Height Hs (m)') -title('Power (kW)') +title('Radiation Capture Width Limit (m)','FontSize',12) colorbar -subplot(1,3,2) +subplot(1,5,4) contourf(T,Hs,p.JPD) xlabel('Wave Period T (s)') ylabel('Wave Height Hs (m)') -title('Probability (%)') +title('Probability at Site (%)','FontSize',12) colorbar -subplot(1,3,3) -contourf(T,Hs,P_matrix .* p.JPD/100) +subplot(1,5,5) +contourf(T,Hs,P_wave .* CW_over_CW_max .* CW_max .* p.JPD) xlabel('Wave Period T (s)') ylabel('Wave Height Hs (m)') -title('Weighted Power (kW)') +title('Weighted Power (W)','FontSize',12) colorbar -improvePlot +%improvePlot end diff --git a/mdocean/plots/sin_desc_fcn_demo.m b/mdocean/plots/sin_desc_fcn_demo.m index c015058f..47330aff 100644 --- a/mdocean/plots/sin_desc_fcn_demo.m +++ b/mdocean/plots/sin_desc_fcn_demo.m @@ -18,7 +18,7 @@ improvePlot set(h,'HandleVisibility','off') -legend('Saturated Signal','Fundamental Amplitude') +legend('Saturated Signal','Fundamental Amplitude','Interpreter','latex') % see what happens when added to phase shifted signal shift_on = false; @@ -49,7 +49,7 @@ drag_actual = sin(x*2*pi) .* abs(sin(x*2*pi)); figure -plot(x,drag_actual,x,drag_DF) +plot(x,drag_actual,'r',x,drag_DF,'b') xlabel('Normalized Time') ylabel('Normalized Drag Force') legend('Squared Signal $\sin(\omega t)|\sin(\omega t)|$',... diff --git a/mdocean/plots/bluewhitered.m b/mdocean/plots/util/bluewhitered.m similarity index 100% rename from mdocean/plots/bluewhitered.m rename to mdocean/plots/util/bluewhitered.m diff --git a/mdocean/plots/hatchfill2.m b/mdocean/plots/util/hatchfill2.m similarity index 100% rename from mdocean/plots/hatchfill2.m rename to mdocean/plots/util/hatchfill2.m diff --git a/mdocean/plots/hatchfill2_demo.m b/mdocean/plots/util/hatchfill2_demo.m similarity index 100% rename from mdocean/plots/hatchfill2_demo.m rename to mdocean/plots/util/hatchfill2_demo.m diff --git a/mdocean/plots/improvePlot.m b/mdocean/plots/util/improvePlot.m similarity index 100% rename from mdocean/plots/improvePlot.m rename to mdocean/plots/util/improvePlot.m diff --git a/mdocean/plots/imrange.m b/mdocean/plots/util/imrange.m similarity index 100% rename from mdocean/plots/imrange.m rename to mdocean/plots/util/imrange.m diff --git a/mdocean/plots/license_bluewhitered.txt b/mdocean/plots/util/license_bluewhitered.txt similarity index 100% rename from mdocean/plots/license_bluewhitered.txt rename to mdocean/plots/util/license_bluewhitered.txt diff --git a/mdocean/plots/license_hatchfill2.txt b/mdocean/plots/util/license_hatchfill2.txt similarity index 100% rename from mdocean/plots/license_hatchfill2.txt rename to mdocean/plots/util/license_hatchfill2.txt diff --git a/mdocean/plots/util/license_table2latex.txt b/mdocean/plots/util/license_table2latex.txt new file mode 100644 index 00000000..8db2cc39 --- /dev/null +++ b/mdocean/plots/util/license_table2latex.txt @@ -0,0 +1,25 @@ +Copyright (c) 2018, Víctor Martínez +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +* Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution +* Neither the name of University of Valladolid nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/mdocean/plots/remove_underscores.m b/mdocean/plots/util/remove_underscores.m similarity index 100% rename from mdocean/plots/remove_underscores.m rename to mdocean/plots/util/remove_underscores.m diff --git a/mdocean/plots/signed_log.m b/mdocean/plots/util/signed_log.m similarity index 100% rename from mdocean/plots/signed_log.m rename to mdocean/plots/util/signed_log.m diff --git a/mdocean/plots/util/table2latex.m b/mdocean/plots/util/table2latex.m new file mode 100644 index 00000000..e56e98aa --- /dev/null +++ b/mdocean/plots/util/table2latex.m @@ -0,0 +1,135 @@ +% ----------------------------------------------------------------------- % +% Function table2latex(T, filename) converts a given MATLAB(R) table into % +% a plain .tex file with LaTeX formatting. % +% % +% Input parameters: % +% - T: MATLAB(R) table. The table should contain only the % +% following data types: numeric, boolean, char or string. +% Avoid including structs or cells. % +% - filename: (Optional) Output path, including the name of the file. +% If not specified, the table will be stored in a % +% './table.tex' file. % +% ----------------------------------------------------------------------- % +% Example of use: % +% LastName = {'Sanchez';'Johnson';'Li';'Diaz';'Brown'}; % +% Age = [38;43;38;40;49]; % +% Smoker = logical([1;0;1;0;1]); % +% Height = [71;69;64;67;64]; % +% Weight = [176;163;131;133;119]; % +% T = table(Age,Smoker,Height,Weight); % +% T.Properties.RowNames = LastName; % +% table2latex(T); % +% ----------------------------------------------------------------------- % +% Version: 1.1 % +% Author: Victor Martinez Cagigal % +% Date: 09/10/2018 % +% E-mail: vicmarcag (at) gmail (dot) com % +% ----------------------------------------------------------------------- % +function table2latex(T, filename, special_col_spec, special_first_row) + + % Error detection and default parameters + if nargin < 2 + filename = 'table.tex'; + fprintf('Output path is not defined. The table will be written in %s.\n', filename); + elseif ~ischar(filename) + error('The output file name must be a string.'); + else + if ~strcmp(filename(end-3:end), '.tex') + filename = [filename '.tex']; + end + end + if nargin < 1, error('Not enough parameters.'); end + if ~istable(T), error('Input must be a table.'); end + if nargin<3 + special_col_spec = []; + end + if nargin<4 + special_first_row = []; + end + + % Parameters + n_col = size(T,2); + col_spec = []; + for c = 1:n_col, col_spec = [col_spec 'l']; end + col_names = strjoin(T.Properties.VariableNames, ' & '); + row_names = T.Properties.RowNames; + if ~isempty(row_names) + col_spec = ['l' col_spec]; + col_names = ['& ' col_names]; + end + + if ~isempty(special_col_spec) + col_spec = strrep(special_col_spec, '\', '\\'); + end + + % Writing header + + table_string = ''; + table_string = append(table_string, ['\\begin{tabular}{' col_spec '}\n']); + + if ~isempty(special_first_row) + table_string = append(table_string, [strrep(special_first_row, '\', '\\') ' \n']); + end + + table_string = append(table_string, [col_names '\\\\ \n']); + table_string = append(table_string, '\\hline \n'); + + + % Writing the data + try + for row = 1:size(T,1) + row_data = cell(1,n_col); + row_data{1,n_col} = []; + for col = 1:n_col + value = T{row,col}; + if isstruct(value) + error('Table must not contain structs.'); + end + while iscell(value) + value = value{1,1}; + end + if ismissing(value) + value = ''; + end + if ~isempty(value) && (ischar(value) || isstring(value)) + value = char(value); + if isstrprop(value(1), 'digit') + value = str2double(value); + end + end + + % Format the output + if isnumeric(value) + if isinf(value) + value = '$\infty$'; + end + exponent = floor(log10(abs(value))/3) * 3; % Round down to nearest multiple of 3 + mantissa = value / 10^exponent; % Calculate mantissa + if exponent==0 + value = sprintf('$%.3g $',mantissa); + else + value = sprintf('$%.3g \\\\cdot 10^{%d}$', mantissa, exponent); + end + end + + row_data{1,col} = char(value); + end + if ~isempty(row_names) + row_data = [row_names{row}, row_data]; + end + row_string = append(strjoin(row_data, ' & '), ' \\\\ \n'); % need 4 \ in row_string and 2 \ in file + table_string = append(table_string, row_string); + clear row_data; + end + catch err + error('Unknown error. Make sure that table only contains chars, strings or numeric values.'); + end + + table_string = append(table_string,'\\end{tabular}'); + + fileID = fopen(filename, 'w'); + fprintf(fileID, table_string); + + % Closing the file + fclose(fileID); +end \ No newline at end of file diff --git a/mdocean/plots/visualize_geometry.m b/mdocean/plots/visualize_geometry.m index 2f0d2059..280828d4 100644 --- a/mdocean/plots/visualize_geometry.m +++ b/mdocean/plots/visualize_geometry.m @@ -61,8 +61,8 @@ function visualize_geometry(x,p,mini,color) set(gcf, 'Color', 'white'); end grid on -ylim([-40 35]) -xlim([-20 20]) +ylim([-40 12]) +xlim([-27 27]) set(waves,'HandleVisibility','off') end @@ -81,5 +81,5 @@ function center_rect(vec,color) function trapezoid(base_1,base_2,y_1,y_2,color) x = [-base_1/2, base_1/2, base_2/2, -base_2/2, -base_1/2]; y = [y_1 y_1 y_2 y_2 y_1]; - plot(x,y,'LineWidth',3,'Color',color) + plot(x,y,'LineWidth',3,'Color',color,'HandleVisibility','off') end diff --git a/mdocean/plots/xdsm.py b/mdocean/plots/xdsm.py index eec2aaea..b8e1edbe 100644 --- a/mdocean/plots/xdsm.py +++ b/mdocean/plots/xdsm.py @@ -6,11 +6,12 @@ x.add_system("opt", OPT, r"\text{Optimizer}") x.add_system("geom", FUNC, r"\text{Geometry}") x.add_system("hydro", FUNC, r"\text{Hydrodynamics}") +x.add_system("solver", SOLVER, r"\text{Iteration}") x.add_system("dynam", FUNC, (r"\text{Dynamics}", r"\text{and Control}")) x.add_system("struct", FUNC, r"\text{Structures}") x.add_system("econ", FUNC, r"\text{Economics}") -x.add_system("F", FUNC, "Objective") -x.add_system("G", FUNC, "Constraints") +x.add_system("F", FUNC, r"\text{Objective}") +x.add_system("G", FUNC, r"\text{Constraints}") x.connect("opt", "geom", (r"\text{Dimensions,}",r"\text{Thicknesses}")) x.connect("opt", "hydro", r"\text{Dimensions}") @@ -18,6 +19,9 @@ x.connect("opt", "struct", (r"\text{Dimensions,}", r"\text{thicknesses}")) x.connect("opt", "econ", (r"\text{Generator}",r"\text{ratings}")) +x.connect("solver","dynam", (r"\text{Dynamic}",r"\text{response guess}")) +x.connect("dynam","solver", (r"\text{Dynamic}",r"\text{response residual}")) + x.connect("geom", "dynam", r"\text{Mass}") x.connect("geom", "econ", r"\text{Material volume}") x.connect("hydro", "dynam", (r"\text{Hydrodynamic}",r"\text{coefficients}")) diff --git a/mdocean/simulation/modules/dynamics/dynamics.m b/mdocean/simulation/modules/dynamics/dynamics.m index 54cbc587..aa4bcde5 100644 --- a/mdocean/simulation/modules/dynamics/dynamics.m +++ b/mdocean/simulation/modules/dynamics/dynamics.m @@ -3,16 +3,17 @@ % use probabilistic sea states for power and PTO force and max amplitude [T,Hs] = meshgrid(in.T,in.Hs); - [P_matrix_mech,X_constraints,B_p,... + [P_matrix_mech,X_constraints_op,B_p,... X_u,X_f,X_s,F_heave_op,... F_surge_op,F_ptrain_max] = get_power_force(in,T,Hs,m_float,m_spar,... - V_d,draft,in.F_max, in.JPD==0); + V_d,draft,in.F_max, ... + in.JPD==0, in.T_f_1); % account for powertrain electrical losses P_matrix_elec = P_matrix_mech * in.eff_pto; % saturate maximum power - P_matrix_elec = min(P_matrix_elec,in.power_max); + P_matrix_elec = min(P_matrix_elec,in.P_max); % weight power across all sea states P_weighted = P_matrix_elec .* in.JPD / 100; @@ -22,10 +23,16 @@ % use max sea states for structural forces % fixme: for storms, should return excitation force, not all hydro force - [~,~,~,~,~,~,F_heave_storm,F_surge_storm,~] = get_power_force(in, ... - in.T_struct, in.Hs_struct, m_float, m_spar, V_d, draft, 0, false(size(in.T_struct))); + [~,X_constraints_storm,~,~,~,~,F_heave_storm,F_surge_storm,~] = get_power_force(in, ... + in.T_struct, in.Hs_struct, m_float, m_spar, V_d, draft, ... + 0, false(size(in.T_struct)), in.T_f_2); F_heave_storm = F_heave_storm * in.F_heave_mult; + % use all X constraints operationally, only use slamming in storm + X_constraints_storm = X_constraints_storm(5:end); + X_constraints_storm = 1 + 0*X_constraints_storm; % fixme this overrides the constraint + X_constraints = [X_constraints_op X_constraints_storm]; + % coefficient of variance (normalized standard deviation) of power P_var = std(P_matrix_elec(:), in.JPD(:)) / P_avg_elec; P_var = P_var * 100; % convert to percentage @@ -34,7 +41,8 @@ end function [P_matrix, X_constraints, B_p, mag_X_u, mag_X_f, mag_X_s,... - F_heave_f, F_surge, F_ptrain_max] = get_power_force(in,T,Hs, m_float, m_spar, V_d, draft, F_max, idx_constraint) + F_heave_f, F_surge, F_ptrain_max] = get_power_force(in,T,Hs, m_float, m_spar, V_d, draft, ... + F_max, idx_constraint, T_f_slam) % get dynamic coefficients for float and spar % fixme: eventually should use in.D_f_in to allow a radial gap between float and spar @@ -48,7 +56,7 @@ in.rho_w, in.g, ... in.use_MEEM, in.harmonics, in.hydro); - X_max = 1e6;%min(Hs / (2*sqrt(2)), in.T_f); + X_u_max = 1e6;%min(Hs / (2*sqrt(2)), in.T_f); % get response: includes drag and force saturation [mag_U,phase_U,... @@ -59,13 +67,14 @@ B_p,K_p] = get_response_drag(w,m_f,m_s,m_c,B_h_f,B_h_s,B_c,K_h_f,K_h_s,... F_f_mag,F_f_phase,F_s_mag,F_s_phase,F_max,... drag_const_f,drag_const_s,mag_v0_f,mag_v0_s, ... - X_max,in.control_type,in.use_multibody,... + X_u_max,in.control_type,in.use_multibody,... in.X_tol,in.phase_X_tol,in.max_drag_iters); % FIXME: check stability of closed loop multibody system - % calculate power - P_matrix = real_P; + % apply empirical fitted "fudge factor" to adjust power from singlebody to multibody + fudge_factor = 23.5 ./ (T.^2 - 14*T + 83) * in.power_scale; + P_matrix = real_P .* fudge_factor; % set values where JPD=0 to 0 to not include them in constraint mag_X_u_const = mag_X_u; @@ -75,27 +84,28 @@ mag_X_s_const = mag_X_s; mag_X_s_const(idx_constraint) = 0; - X_max = max(mag_X_u_const,[],'all'); + X_u_max = max(mag_X_u_const,[],'all'); + X_f_max = max(mag_X_f_const,[],'all'); + % extra height on spar after accommodating float displacement - h_s_extra_up = (in.h_s - in.T_s - (in.h_f - in.T_f_2) - X_max) / in.h_s; - h_s_extra_down = (in.T_s - in.T_f_2 - X_max) / in.h_s; + h_s_extra_up = (in.h_s - in.T_s - (in.h_f - in.T_f_2) - X_u_max) / in.h_s; + h_s_extra_down = (in.T_s - in.T_f_2 - X_u_max) / in.h_s; % sufficient length of float support tube - h_fs_extra = in.h_fs_clear / X_max - 1; + h_fs_extra = in.h_fs_clear / X_u_max - 1; % prevent violation of linear wave theory X_max_linear = 1/10 * in.D_f; - X_below_linear = X_max_linear / X_max - 1; + X_below_linear = X_max_linear / X_f_max - 1; % prevent rising out of the water wave_amp = Hs/(2*sqrt(2)); - X_over_A = mag_X_u_const ./ wave_amp; - T_over_A = in.T_f_2 ./ wave_amp; - R = T_over_A ./ sqrt(1 + X_over_A.^2 - 2 * X_over_A .* cos(phase_X_u)); % ratio derived on p88-89 of notebook - long_draft = R-1; - small_diameter = 2*pi - 2*real(acos(R)) - k_wvn * in.D_f; - X_below_wave = max(long_draft, small_diameter); % one or the other is required, but not necessarily both + theta_slam = max(0, -k_wvn * in.D_f / 2 + abs(pi - phase_X_f)); + X_slam = sqrt( T_f_slam^2 - (wave_amp .* sin(theta_slam)).^2 ) - wave_amp .* cos(theta_slam); + X_slam( imag(X_slam)~=0 ) = 0; % slamming occurs even for stationary body + X_below_wave = X_slam ./ mag_X_u_const - 1; + X_below_wave(~isfinite(X_below_wave)) = 1; % constraint always satisfied when JPD=0 X_constraints = [h_s_extra_up, h_s_extra_down, h_fs_extra, X_below_linear, X_below_wave(:).']; diff --git a/mdocean/simulation/modules/econ.m b/mdocean/simulation/modules/econ.m deleted file mode 100644 index 2a089505..00000000 --- a/mdocean/simulation/modules/econ.m +++ /dev/null @@ -1,34 +0,0 @@ - -function [LCOE, capex_design_dependent_perwec, ... - capex, opex] = econ(m_m, M, cost_m, N_WEC, P_elec, FCR, ... - cost_perN, cost_perW, F_max, P_max, efficiency) - -structural_cost = m_m.* cost_m; -devicestructure = structural_cost(M); - -% costs taken from 'CBS (Total)' tab of the RM3 cost breakdown structure -% https://catalog.data.gov/ne/dataset/reference-model-3-cost-breakdown-rm3-wave-point-absorber -development = 4553000; -infrastructure = 990000; -mooring = N_WEC * 525000; -pto = 623000 + F_max * cost_perN + P_max * cost_perW; -profitmargin = 356000; -installation = 5909000; -contingency = 1590000; - -capex_design_dependent_perwec = devicestructure + pto; -capex = development + infrastructure + mooring + capex_design_dependent_perwec * N_WEC ... - + profitmargin + installation + contingency; - -% opex reflects cost of operation, postinstall, replacement, consumables, and -% insurance -% opex equation was generated by Casio online power regression calculator -opex = 1193000 * N_WEC^0.4433; - -hr_per_yr = 8766; -P_avg = N_WEC * P_elec * efficiency; -aep = P_avg * hr_per_yr / 1000; % annual energy production: W to kWh per year - -LCOE = (FCR*capex + opex)/aep; - -end \ No newline at end of file diff --git a/mdocean/simulation/modules/econ/LCOE_from_capex_design_power.m b/mdocean/simulation/modules/econ/LCOE_from_capex_design_power.m new file mode 100644 index 00000000..23a7c6f7 --- /dev/null +++ b/mdocean/simulation/modules/econ/LCOE_from_capex_design_power.m @@ -0,0 +1,28 @@ +function [LCOE,capex,opex] = LCOE_from_capex_design_power(capex_design_dep, N_WEC, P_elec, FCR, efficiency) + + % costs taken from 'CBS (Total)' tab of the RM3 cost breakdown structure + % https://catalog.data.gov/ne/dataset/reference-model-3-cost-breakdown-rm3-wave-point-absorber + + % design-independent cost per wec + % includes development, infrastructure, mooring, profitmargin, + % installation, contingency. (no decommissioning intentionally) + alpha_non_design = 0.741; + capex_non_design_dep = 12.68e6 * N_WEC^(-alpha_non_design) + 1.24e6; + + % total capex + capex_per_wec = capex_design_dep + capex_non_design_dep; + capex = capex_per_wec * N_WEC; + + % opex = operation, postinstall, replacement, consumables, and insurance + alpha_opex = 0.5567; + opex_per_wec = 1.193e6 * N_WEC^-alpha_opex; + opex = opex_per_wec * N_WEC; + + % power and energy + hr_per_yr = 8766; + P_avg = N_WEC * P_elec * efficiency; + AEP = P_avg * hr_per_yr / 1000; % annual energy production: W to kWh per year + + LCOE = (FCR*capex + opex)./AEP; + +end \ No newline at end of file diff --git a/mdocean/simulation/modules/econ/econ.m b/mdocean/simulation/modules/econ/econ.m new file mode 100644 index 00000000..b1eab66d --- /dev/null +++ b/mdocean/simulation/modules/econ/econ.m @@ -0,0 +1,35 @@ + +function [LCOE, capex_design_dep, ... + capex, opex, pto, devicestructure] = econ(mass_material, M, cost_perkg_mult, N_WEC, P_elec, FCR, ... + cost_perN_mult, cost_perW_mult, F_max, P_max, efficiency) + + [capex_design_dep, ... + pto, devicestructure] = design_cost_model(mass_material, M, cost_perkg_mult, N_WEC, ... + cost_perN_mult, cost_perW_mult, F_max, P_max); + + [LCOE,capex,opex] = LCOE_from_capex_design_power(capex_design_dep, N_WEC, P_elec, FCR, efficiency); + +end + +function [capex_design_dep, ... + pto, devicestructure] = design_cost_model(mass_material, M, cost_perkg_mult, N_WEC, ... + cost_perN_mult, cost_perW_mult, F_max, P_max) + + % from RM3 CBS 1.4 and 1.5, with curve fits done in dev/design_cost_scaling.m + + % structural cost per wec + alpha_struct = 0.481; + cost_per_kg = ( 1.64e6 + 1.31e6 * N_WEC^(-alpha_struct) ) / 687000 * cost_perkg_mult(M); + devicestructure = cost_per_kg * mass_material; + + % PTO cost per wec + alpha_pto = 0.206; + pto_const = 92593 + 1051 *N_WEC^(-alpha_pto); + pto_power = ( 0.4454 + 0.9099*N_WEC^(-alpha_pto) ) * P_max * cost_perN_mult; + pto_force = ( 0.0648 + 0.0893*N_WEC^(-alpha_pto) ) * F_max * cost_perW_mult; + pto = pto_const + pto_power + pto_force; + + % sum design-dependent cost per wec + capex_design_dep = devicestructure + pto; + +end \ No newline at end of file diff --git a/mdocean/simulation/modules/geometry.m b/mdocean/simulation/modules/geometry.m index 104d975d..c7c43c93 100644 --- a/mdocean/simulation/modules/geometry.m +++ b/mdocean/simulation/modules/geometry.m @@ -1,11 +1,12 @@ function [V_d, m_m, m_f_tot, m_s_tot,... A_c, A_lat_sub, ... I, T, V_f_pct, V_s_pct, GM,... - A_dt, L_dt, t_d, mass,... + A_dt, L_dt, mass,... CB_f_from_waterline,CG_f_from_waterline] = geometry(D_s, D_f, D_f_in, D_f_b, ... T_f_1, T_f_2, h_f, h_s, h_fs_clear, D_f_tu, ... t_f_t, t_f_r, t_f_c, t_f_b, t_f_tu, t_s_r, t_d_tu, ... - D_d, D_d_tu, theta_d_tu, T_s, h_d, t_d_max, ... + D_d, D_d_tu, theta_d_tu, T_s, h_d, t_d, ... + h_stiff_f, w_stiff_f, num_sect_f, ... h_stiff_d, w_stiff_d, num_stiff_d, ... M, rho_m, rho_w, m_scale) @@ -52,7 +53,7 @@ % t_f_tu - radial thicknss of float support tube walls %% Float -num_gussets = 24; +num_gussets = 2*num_sect_f; num_gussets_loaded_lateral = 2; % float cross sectional and lateral area for structural purposes @@ -78,8 +79,14 @@ L_ft = sqrt(horiz_leg_tube^2 + vert_leg_tube^2); % length of float tube V_f_tubes = num_float_tubes * (D_f_tu^2 - (D_f_tu - 2*t_f_tu)^2) * pi/4 * L_ft; -V_sf_m = V_top_plate + V_bot_plate + V_outer_rim ... - + V_inner_rim + V_bot_slant + V_gussets + V_f_tubes; +% float stiffeners +num_stiff_per_sect_f = 2; % top and bottom - neglecting the side stiffeners +A_stiff_f = num_stiff_per_sect_f * num_sect_f * h_stiff_f * w_stiff_f; +len_stiff_f = (D_f - D_f_in)/2 * 0.75; % goes ~3/4 of the way along the float +V_stiff_f = len_stiff_f * A_stiff_f; + +V_sf_m = V_top_plate + V_bot_plate + V_outer_rim + V_inner_rim ... + + V_bot_slant + V_gussets + V_f_tubes + V_stiff_f; m_f_m = V_sf_m * rho_m(M) * m_scale; % mass of float material without ballast @@ -111,7 +118,15 @@ % vertical column material use D_vc_i = D_s - 2 * t_s_r; % spar column inner diameter A_vc_c = pi/4 * (D_s^2 - D_vc_i^2); % spar column cross sectional area -V_vc_m = A_vc_c * (h_s - h_d); % volume of column material +V_vc_m_body = A_vc_c * (h_s - h_d); +A_vc_caps = pi/4 * D_vc_i^2; +t_vc_caps = (1/2 + 2.5) * 0.0254; % middle is 2.5", top is 0.5" +V_vc_caps = A_vc_caps * t_vc_caps; +num_stiff_vc = 12; +A_vc_stiff = 0.658 + 2*0.652 + 0.350; % triangles +t_vc_stiff = 0.0254; +V_vc_m_stiff = num_stiff_vc * A_vc_stiff * t_vc_stiff; +V_vc_m = V_vc_m_body + V_vc_m_stiff + V_vc_caps; % volume of column material % damping plate material use A_d = pi/4 * D_d^2; % damping plate itself @@ -119,7 +134,6 @@ L_dt = (D_d - D_s) / (2*cos(theta_d_tu)); D_dt_i = D_d_tu - 2 * t_d_tu; A_dt = pi/4 * (D_d_tu^2 - D_dt_i^2); % support tube area -t_d = min(h_d, t_d_max); % material thickenss of damping plate num_unique_stiffeners = length(h_stiff_d)/2; num_stiff_repeats = num_stiff_d / num_unique_stiffeners; A_stiff_d = num_stiff_repeats * sum(h_stiff_d .* w_stiff_d); diff --git a/mdocean/simulation/run/check_max_CW.m b/mdocean/simulation/run/check_max_CW.m index c06b51e5..152c1453 100644 --- a/mdocean/simulation/run/check_max_CW.m +++ b/mdocean/simulation/run/check_max_CW.m @@ -1,17 +1,10 @@ -function [ratio] = check_max_CW(filename_uuid) +function [ratio, P_wave, CW_max] = check_max_CW(filename_uuid) p = parameters(); - p.cost_m = [0 0 0]; % hack that makes cost constant, so minimizing LCOE is actually maximizing power - p.power_max = Inf; - b = var_bounds(); - b.filename_uuid = filename_uuid; - x0 = b.X_start_struct; + %b.filename_uuid = filename_uuid; - % run LCOE minimization (effectively power maximization due to hack above) - X_opt = gradient_optim(x0,p,b,1); + [~,val] = max_avg_power(p,b); - % plug back into simulation to get unsaturated power - [~, ~, ~, ~, val] = simulation(X_opt,p); P_mech = val.P_mech; % calculate capture width diff --git a/mdocean/simulation/run/is_feasible.m b/mdocean/simulation/run/is_feasible.m index b14c931f..e151e2b8 100644 --- a/mdocean/simulation/run/is_feasible.m +++ b/mdocean/simulation/run/is_feasible.m @@ -1,7 +1,15 @@ -function [feasible, failed] = is_feasible(g_nonlin, x, p, b) +function [feasible, failed, feasible_lin] = is_feasible(g_nonlin, x, p, b, idx_ignore) + +if nargin<5 + idx_ignore = false(size(g_nonlin)); +end tol = -.01; -feasible_nonlin = all(g_nonlin >= tol); +feasible_nonlin = all(g_nonlin(~idx_ignore) >= tol); + +% The linear constraint currently assume a certain design variable order. +% If you change the order, this assert reminds you to update A_ineq & b_ineq. +assert(all( strcmp(b.var_names(1:4),{'D_s','D_f','T_f_2','h_s'}) )) A_ineq = [-p.D_d_over_D_s 0 0 0; 1 -1 0 0; diff --git a/mdocean/simulation/run/validate_nominal_RM3.m b/mdocean/simulation/run/validate_nominal_RM3.m index f3f9a17e..0dd3f858 100644 --- a/mdocean/simulation/run/validate_nominal_RM3.m +++ b/mdocean/simulation/run/validate_nominal_RM3.m @@ -1,8 +1,8 @@ function [feasible,failed,simulated,actual,tab,fig] = validate_nominal_RM3(mode) p = parameters(mode); p.N_WEC = 1; - p.power_max = 286000; p.LCOE_max = 10; % set large max LCOE to avoid failing feasibility check + p.control_type = 'damping'; b = var_bounds(mode); X = [b.X_noms; 1]; @@ -10,13 +10,15 @@ [~, ~, ~, g, simulated] = simulation(X,p); % whether nominal violates constraints - [feasible,failed] = is_feasible(g, X, p, b); + idx_ignore = strcmp(b.constraint_names,'irrelevant_max_force'); + [feasible,failed] = is_feasible(g, X, p, b, idx_ignore); % comparison of simulated and actual values if nargout > 2 actual = validation_inputs(mode); fig = figure; - t = tiledlayout(fig,1,3); + econ_fields = {'capex','opex','LCOE','capex_design','capex_struct','capex_PTO','J_capex_design'}; + t = tiledlayout(fig,1,length(econ_fields)); fields = fieldnames(actual); % for economic validation, sweep N_WEC @@ -32,7 +34,7 @@ field = fields{i}; % if the field is economic, plot vs N_WEC - if any(strcmp(field,{'capex','opex','LCOE'})) + if any(strcmp(field,econ_fields)) simulated.(field) = [simulated_diff_N_WEC.(field)]; ax = nexttile(t); diff --git a/mdocean/simulation/simulation.m b/mdocean/simulation/simulation.m index 2b5c6030..cecfffa0 100644 --- a/mdocean/simulation/simulation.m +++ b/mdocean/simulation/simulation.m @@ -5,21 +5,29 @@ %% Assemble inputs in = p; -in.D_s = X(1); % inner diameter of float (m) -in.D_f = X(2); % outer diameter of float (m) -in.T_f_2 = X(3); % draft of float (m) -in.h_s = X(4); % total height of spar (m) -in.F_max = X(5)*1e6; % max powertrain force (N) -in.B_p = X(6)*1e6; % controller (powertrain) damping (Ns/m) -in.h_fs_clear = X(7); % vertical clearance between float tubes and spar when at rest (m) -in.t_f_t = X(8)*1e-3; % float top thickness (m) -in.t_f_r = X(9)*1e-3; % float radial wall thickness (m) -in.t_f_c = X(10)*1e-3; % float circumferential gusset thickness (m) -in.t_f_b = X(11)*1e-3; % float bottom thickness (m) -in.t_s_r = X(12)*1e-3; % vertical column thickness (m) -in.t_d_tu = X(13)*1e-3; % damping plate support tube radial wall thickness (m) -in.P_max = X(14)*1e3; % maximum power (W) -in.M = X(15); % material (-) +in.D_s = X(1); % inner diameter of float (m) +in.D_f = X(2); % outer diameter of float (m) +in.T_f_2 = X(3); % draft of float (m) +in.h_s = X(4); % total height of spar (m) +in.h_fs_clear = X(5); % vertical clearance between float tubes and spar when at rest (m) +in.F_max = X(6) * 1e6; % max powertrain force (N) +in.P_max = X(7) * 1e3; % maximum power (W) +in.t_f_b = X(8) * 1e-3; % float bottom thickness (m) +in.t_s_r = X(9) * 1e-3; % vertical column thickness (m) +in.t_d = X(10)* 1e-3; % damping plate support tube radial wall thickness (m) +in.h_stiff_f = X(11); % float stiffener height (m) +in.h1_stiff_d = X(12); % damping plate stiffener larger height (m) +in.M = X(13); % material index (-) + +% float thickness ratios +in.t_f_r = in.t_f_b * p.t_f_r_over_t_f_b; % float radial wall thickness (m) +in.t_f_c = in.t_f_b * p.t_f_c_over_t_f_b; % float circumferential gusset thickness (m) +in.t_f_t = in.t_f_b * p.t_f_t_over_t_f_b; % float bottom thickness (m) + +% stiffener thickness ratios +in.w_stiff_f = p.w_over_h_stiff_f * in.h_stiff_f; +in.h_stiff_d = p.h_over_h1_stiff_d * in.h1_stiff_d; +in.w_stiff_d = p.w_over_h1_stiff_d * in.h1_stiff_d; % Geometric similarity to maintain constant damping ratio % D_s sets D_d, T_s, h_d @@ -42,14 +50,15 @@ [V_d, m_m, m_f_tot, m_s_tot, ... A_c, A_lat_sub, ... I, T, V_f_pct, V_s_pct, GM, ... - A_dt, L_dt, t_d] = geometry(in.D_s, in.D_f, in.D_f_in, in.D_f_b, ... - in.T_f_1, in.T_f_2, in.h_f, in.h_s, ... - in.h_fs_clear, in.D_f_tu, in.t_f_t, ... - in.t_f_r, in.t_f_c, in.t_f_b, in.t_f_tu, ... - in.t_s_r, in.t_d_tu, in.D_d, in.D_d_tu, ... - in.theta_d_tu, in.T_s, in.h_d, in.t_d_max,... - in.h_stiff_d, in.w_stiff_d, in.num_stiff_d, ... - in.M, in.rho_m, in.rho_w, in.m_scale); + A_dt, L_dt] = geometry(in.D_s, in.D_f, in.D_f_in, in.D_f_b, ... + in.T_f_1, in.T_f_2, in.h_f, in.h_s, ... + in.h_fs_clear, in.D_f_tu, in.t_f_t, ... + in.t_f_r, in.t_f_c, in.t_f_b, in.t_f_tu, ... + in.t_s_r, in.t_d_tu, in.D_d, in.D_d_tu, ... + in.theta_d_tu, in.T_s, in.h_d, in.t_d,... + in.h_stiff_f, in.w_stiff_f, in.num_sections_f, ... + in.h_stiff_d, in.w_stiff_d, in.num_stiff_d, ... + in.M, in.rho_m, in.rho_w, in.m_scale); m_f_tot = max(m_f_tot,1e-3); % zero out negative mass produced by infeasible inputs @@ -61,17 +70,20 @@ [FOS_float,FOS_spar,FOS_damping_plate,... FOS_spar_local] = structures(... F_heave_storm, F_surge_storm, F_heave_op, F_surge_op, ... % forces - in.h_s, in.T_s, in.D_s, in.D_f, in.D_f_in, in.num_sections, in.D_f_tu, in.D_d, L_dt, in.theta_d_tu,in.D_d_tu,... % bulk dimensions - in.t_s_r, I, A_c, A_lat_sub, in.t_f_b, in.t_f_t, t_d, in.t_d_tu, in.h_d, A_dt, ... % structural dimensions + in.h_s, in.T_s, in.D_s, in.D_f, in.D_f_in, in.num_sections_f, ... % bulk dimensions + in.D_f_tu, in.D_d, L_dt, in.theta_d_tu,in.D_d_tu,... % more bulk dimensions + in.t_s_r, I, A_c, A_lat_sub, in.t_f_b, in.t_f_t, in.t_d, in.t_d_tu, in.h_d, A_dt, ... % structural dimensions in.h_stiff_f, in.w_stiff_f, in.h_stiff_d, in.w_stiff_d,... % stiffener thicknesses - in.M, in.rho_w, in.g, in.sigma_y, in.sigma_e, in.E, in.nu, ... + in.M, in.rho_w, in.g, in.sigma_y, in.sigma_e, in.E, in.nu, ... % constants in.num_terms_plate, in.radial_mesh_plate, in.num_stiff_d); -[LCOE,capex_design] = econ(m_m, in.M, in.cost_m, in.N_WEC, P_avg_elec, in.FCR, in.cost_perN, in.cost_perW, in.F_max, in.P_max, in.eff_array); +[LCOE,capex_design] = econ(m_m, in.M, in.cost_perkg_mult, in.N_WEC, P_avg_elec, ... + in.FCR, in.cost_perN_mult, in.cost_perW_mult, ... + in.F_max, in.P_max, in.eff_array); J_capex_design = capex_design / 1e6; % convert $ to $M %% Assemble constraints g(x) >= 0 -num_g = 23+numel(p.JPD); +num_g = 20+numel(p.JPD)+length(p.T_struct); g = zeros(1,num_g); g(1) = V_f_pct; % prevent float too heavy g(2) = 1 - V_f_pct; % prevent float too light @@ -80,27 +92,24 @@ g(5) = GM; % pitch stability of float-spar system g(6) = FOS_float(1) / p.FOS_min - 1; % float survives max force g(7) = FOS_float(2) / p.FOS_min - 1; % float survives fatigue -g(8) = in.t_f_t / in.t_f_b - 0.5 / 0.56; % float top thickness ratio placeholder -g(9) = in.t_f_r / in.t_f_b - 0.44 / 0.56; % float radial thickness ratio placeholder -g(10) = in.t_f_c / in.t_f_b - 0.44 / 0.56; % float circumferential thickness ratio placeholder -g(11) = FOS_spar(1) / p.FOS_min - 1; % spar survives max force -g(12) = FOS_spar(2) / p.FOS_min - 1; % spar survives fatigue -g(13) = FOS_damping_plate(1) * in.FOS_mult_d / p.FOS_min - 1; % damping plate survives max force -g(14) = FOS_damping_plate(2) * in.FOS_mult_d / p.FOS_min - 1; % damping plate survives fatigue -g(15) = FOS_spar_local(1) / p.FOS_min - 1; % spar survives max force in local buckling -g(16) = FOS_spar_local(2) / p.FOS_min - 1; % spar survives fatigue in local buckling -g(17) = P_avg_elec; % positive power +g(8) = FOS_spar(1) / p.FOS_min - 1; % spar survives max force +g(9) = FOS_spar(2) / p.FOS_min - 1; % spar survives fatigue +g(10) = FOS_damping_plate(1) * in.FOS_mult_d / p.FOS_min - 1; % damping plate survives max force +g(11) = FOS_damping_plate(2) * in.FOS_mult_d / p.FOS_min - 1; % damping plate survives fatigue +g(12) = FOS_spar_local(1) / p.FOS_min - 1; % spar survives max force in local buckling +g(13) = FOS_spar_local(2) / p.FOS_min - 1; % spar survives fatigue in local buckling +g(14) = P_avg_elec; % positive power %1 + min(Kp_over_Ks,[],'all'); % spar heave stability (positive effective stiffness) -g(18) = p.LCOE_max/LCOE - 1; % prevent more expensive than threshold +g(15) = p.LCOE_max/LCOE - 1; % prevent more expensive than threshold %g(19) = P_avg_elec/p.avg_power_min - 1; % prevent less avg power than threshold -g(19) = F_ptrain_max/in.F_max - 1; % prevent irrelevant max force - +g(16) = F_ptrain_max/in.F_max - 1; % prevent irrelevant max force - % this constraint should always be active % and is only required when p.cost_perN = 0. -g(20) = X_constraints(1); % prevent float rising above top of spar -g(21) = X_constraints(2); % prevent float going below bottom of spar -g(22) = X_constraints(3); % prevent float support tube (PTO attachment) from hitting spar -g(23) = X_constraints(4); % float amplitude obeys linear theory -g(24:end) = X_constraints(5:end); % prevent rising out of water/slamming +g(17) = X_constraints(1); % prevent float rising above top of spar +g(18) = X_constraints(2); % prevent float going below bottom of spar +g(19) = X_constraints(3); % prevent float support tube (PTO attachment) from hitting spar +g(20) = X_constraints(4); % float amplitude obeys linear theory +g(21:end) = X_constraints(5:end); % prevent rising out of water/slamming criteria = all(~isinf([g LCOE P_var])) && all(~isnan([g LCOE P_var])) && all(isreal([g LCOE P_var])); if ~criteria @@ -108,23 +117,28 @@ end if nargout > 4 % if returning extra struct output for validation - [~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~,... + [~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~,... mass,CB_f,CG_f] = geometry(in.D_s, in.D_f, in.D_f_in, in.D_f_b, ... in.T_f_1, in.T_f_2, in.h_f, in.h_s, ... in.h_fs_clear, in.D_f_tu, in.t_f_t, ... in.t_f_r, in.t_f_c, in.t_f_b, in.t_f_tu, ... in.t_s_r, in.t_d_tu, in.D_d, in.D_d_tu, ... - in.theta_d_tu, in.T_s, in.h_d, in.t_d_max,... + in.theta_d_tu, in.T_s, in.h_d, in.t_d,... + in.h_stiff_f, in.w_stiff_f, in.num_sections_f, ... in.h_stiff_d, in.w_stiff_d, in.num_stiff_d, ... in.M, in.rho_m, in.rho_w, in.m_scale); - [~,~,capex,opex] = econ(m_m, in.M, in.cost_m, in.N_WEC, P_avg_elec, in.FCR, ... - in.cost_perN, in.cost_perW, in.F_max, in.P_max, in.eff_array); - [~, ~, ~, ~, ~, ~, ~, B_p,X_u,X_f,X_s,P_matrix_mech] = dynamics(in, m_f_tot, m_s_tot, V_d, T); + [~,~,capex,opex,pto, devicestructure] = econ(m_m, in.M, in.cost_perkg_mult, in.N_WEC, P_avg_elec, in.FCR, ... + in.cost_perN_mult, in.cost_perW_mult, in.F_max, in.P_max, in.eff_array); + [~, ~, ~, ~, ~, ~, ~, ~, ~, B_p,X_u,X_f,X_s,P_matrix_mech] = dynamics(in, m_f_tot, m_s_tot, V_d, T); val.mass_f = mass(1); val.mass_vc = mass(2); val.mass_rp = mass(3); val.mass_tot = sum(mass); val.capex = capex; + val.capex_design = capex_design; + val.J_capex_design = J_capex_design; + val.capex_struct = devicestructure; + val.capex_PTO = pto; val.opex = opex; val.LCOE = LCOE; val.power_avg = P_avg_elec; diff --git a/tests/test.m b/tests/test.m index a5bb8fe3..a1e798ec 100644 --- a/tests/test.m +++ b/tests/test.m @@ -5,7 +5,7 @@ properties (Constant) run_slow_tests = false; - slow_figs = [16:21 24 25]; + slow_figs = [16:21]; slow_tabs = 7; end @@ -29,8 +29,8 @@ properties (TestParameter) field_report = fieldnames(validation_inputs('report')); field_wecsim = fieldnames(validation_inputs('wecsim')); - rel_tol_report = {.1,.1,.1,.1,.01,.01,.25,.25,.25,.1,.1,.1,.1,.1}; - rel_tol_wecsim = {.01,.01,.01,.01, 0.1,0.1}; + rel_tol_report = {.1,.1,.1,.1,.01,.01,.25,.25,.25,.1,.1,.1,.1,.1,.1,.1,.1,.1}; + rel_tol_wecsim = {.01,.01,.01,.01, 0.1,0.1,.1,.1}; which_figs = test.enumerateFigs() which_tabs = test.enumerateTabs() end @@ -132,7 +132,7 @@ function deleteGeneratedFiles(testCase) % run every figure and log it function allFiguresRun(testCase, which_figs, which_tabs) - [success_criterion,fig_out,tab_out] = all_figures(which_figs,which_tabs,testCase.uuid.Value); + [success_criterion,fig_out,tab_out, num_figs] = all_figures(which_figs,which_tabs,testCase.uuid.Value); if isempty(success_criterion) success_criterion = 1;