diff --git a/mdocean/optimization/gradient_optim.m b/mdocean/optimization/gradient_optim.m index a24bc4b4..207a3fdf 100644 --- a/mdocean/optimization/gradient_optim.m +++ b/mdocean/optimization/gradient_optim.m @@ -1,10 +1,10 @@ -function [Xs_opt, objs_opt, flags, probs] = gradient_optim(x0_input,p,b,which_objs) +function [Xs_opt, objs_opt, flags, probs, lambda, gs] = gradient_optim(x0_input,p,b,which_objs) if nargin == 0 % set default parameters if function is run without input clc;close all p = parameters(); - b = var_bounds(); + b = var_bounds(p); x0_input = b.X_start_struct; display = 'iter'; plotfn = @optimplotfval; @@ -40,14 +40,14 @@ for matl = 1%1:2:3 %b.M_min : b.M_max X = [D_f D_s_ratio h_f_ratio T_s_ratio F_max B_p w_n matl]; - [Xs_opt, objs_opt, flags, probs] = optimize_both_objectives(X,p,b,x0_input,opts,ploton,which_objs); + [Xs_opt, objs_opt, flags, probs,lambda,g,gs] = 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,lambda,g,gs] = optimize_both_objectives(X,p,b,x0_input,opts,ploton,which_objs) [LCOE, P_var, ~, g] = fcn2optimexpr(@simulation,X,p,... 'OutputSize',{[1,1],[1,1],size(p.JPD),[1, 14]},... @@ -98,9 +98,10 @@ end X_opt = [X_opt_raw; evaluate(X(8),struct())]; % add material back onto design vector - [out(1),out(2)] = simulation(X_opt,p); % rerun sim + [out(1),out(2),~,g] = simulation(X_opt,p); % rerun sim assert(out(which_obj) == obj_opt) % check correct reordering of X_opt elements + gs(:,i) = g; Xs_opt(:,i) = X_opt; objs_opt(i) = obj_opt; flags(i) = flag; @@ -118,5 +119,4 @@ array2table(table_data,'RowNames',b.var_names(1:end-1),... 'VariableNames',{'Min LCOE','Min cv','Min bound','Max bound'}) end - -end +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 7c02558c..97e9bee8 100644 --- a/mdocean/optimization/multiobjective/pareto_curve_heuristics.m +++ b/mdocean/optimization/multiobjective/pareto_curve_heuristics.m @@ -9,6 +9,9 @@ function pareto_curve_heuristics() X = [X ones(length(X),1)]; % add eigth column for material LCOE = fval(:,1); Pvar = fval(:,2); + lambdaActive = lambda.active'; + lambdaLower = lambda.lower'; + lambdaUpper = lambda.upper'; [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,... @@ -39,6 +42,8 @@ function pareto_curve_heuristics() design_heuristics_plot(LCOE, minLCOE, idx_best_LCOE, x_best_LCOE, ... Pvar, minPvar, idx_best_Pvar, X, idxo, p.LCOE_max) + %% plots Lagrange multipliers as a fn of percent along the pareto + lagrange_multiplier_plot(lambdaActive,lambdaUpper,lambdaLower) end %% function [LCOE, minLCOE, idx_best_LCOE, LCOE_nom, ... @@ -70,7 +75,7 @@ function pareto_curve_heuristics() P_var_solar = 125; % balanced design - [~,idx_balanced] = min(abs(Pvar-100)); + [~,idx_balanced] = min(abs(Pvar-40)); LCOE_balanced = LCOE(idx_balanced); P_var_balanced = Pvar(idx_balanced); @@ -126,32 +131,32 @@ function pareto_curve_heuristics() 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(LCOE_balanced+.02,P_var_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]) + axes('Position',[.28 .7 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]) + axes('Position',[.52 .15 mini_plot_size]) box on visualize_geometry(x_best_Pvar,p,true); set(gca,'XTickLabel',[],'YTickLabel',[]) % balanced - axes('Position',[.10 .28 mini_plot_size]) + axes('Position',[.22 .25 mini_plot_size]) box on visualize_geometry(x_balanced,p,true); set(gca,'XTickLabel',[],'YTickLabel',[]) % RM3 - axes('Position',[.7 .53 mini_plot_size]) + axes('Position',[.7 .5 mini_plot_size]) box on visualize_geometry(x_nom,p,true); set(gca,'XTickLabel',[],'YTickLabel',[]) @@ -253,4 +258,78 @@ function pareto_curve_heuristics() improvePlot cent = char(0162); legend(['LCOE (' cent '/kWh)'],'c_v (%)',Location='northeast') +end + +%% +function [] = lagrange_multiplier_plot(lambdaActive, ... + lambdaUpper,lambdaLower) + figure + hold on + [m,n]=size(lambdaActive); + color_cell = {'.r','.g','.b','.c','.m','.y','.k',"#FF0000",'ob','+b','*b','xb','-b','|b'}; + %color_cell = {'r','g','b','c','m','y','k',[0 0.4470 0.7410],[0.8500 0.3250 0.0980],[0.9290 0.6940 0.1250],[0.4940 0.1840 0.5560],[0.4660 0.6740 0.1880],[0.3010 0.7450 0.9330],[0.6350 0.0780 0.1840]}; + for i=1:m + for j = 1:n + %color = color_cell{rem(j,15)}; + plot(i*100/60,lambdaActive(i,j),'.','MarkerSize',12); + end + end + %plot(pct_angle,lambda_active_sorted) + title("Lagrange Multipliers",'FontSize',20) + xlabel("Percent Along the Pareto Front",'FontSize',15) + ylabel("Lagrange Multiplier",'FontSize',15) + legend("prevent float too heavy","prevent float too light", ... + "prevent spar too heavy","prevent spar too light","stability",... + "float survives max force","spar survives max force",... + "damping plate survives max force","spar doesn't buckle", ... + 'positive power','damping plate diameter','prevent float rising above top of spar',... + 'prevent too expensive','prevent irrelevant max force') + hold off + + % lower bound plot + figure + hold on + [m,n]=size(lambdaLower); + % I'll change these later + for i=1:m + for j = 1:n + color = color_cell{rem(j,15)}; + plot(i*100/60,lambdaLower(i,j),color,'MarkerSize',12); + end + end + title('Lower Bound Active Lagrange Multipliers') + xlabel('Percent Along the Pareto Curve') + ylabel('Lagrange Multiplier') + xlim([-1,101]); + ylim([-.05,.2]) + legend('WEC Surface Float Outer Diameter', ... + 'Ratio of WEC Surface Float Inner Diameter to Outer Diameter', ... + 'Ratio of WEC Surface Float Height to Outer Diameter', ... + 'Percent of WEC Spar Submergence','Maximum Powertrain Force', ... + 'Powertrain/Controller Damping','Controller Natural Frequency') + hold off + + % upper bound plot + figure + hold on + [m,n]=size(lambdaUpper) + % I'll change these later + for i=1:m + for j = 1:n + color = color_cell{rem(j,15)}; + plot(i*100/60,lambdaUpper(i,j),color,'MarkerSize',12); + end + end + title('Upper Bound Active Lagrange Multipliers') + xlabel('Percent Along the Pareto Curve') + ylabel('Lagrange Multiplier') + %xlim([-1,61]); + ylim([-.02,.04]) + legend('WEC Surface Float Outer Diameter', ... + 'Ratio of WEC Surface Float Inner Diameter to Outer Diameter', ... + 'Ratio of WEC Surface Float Height to Outer Diameter', ... + 'Percent of WEC Spar Submergence','Maximum Powertrain Force', ... + 'Powertrain/Controller Damping','Controller Natural Frequency') + improvePlot + hold off end \ No newline at end of file diff --git a/mdocean/optimization/multiobjective/pareto_search.m b/mdocean/optimization/multiobjective/pareto_search.m index ab5bbd50..a8445d1d 100644 --- a/mdocean/optimization/multiobjective/pareto_search.m +++ b/mdocean/optimization/multiobjective/pareto_search.m @@ -60,10 +60,41 @@ idx = constraint_active_plot(residuals,fval,tol); cols = [1 3 6 5 4 2 7]; - x_sorted = x(idx,cols) + x_sorted = x(idx,cols); + + % solve for lambda values + [rows,~] = size(x); + for i = 1:length(LCOE_seeds) + p.LCOE_max = LCOE_seeds(i); + which_objs = 2; + end + + for i = 1:rows + input = x(i,:); + x_0_new.D_f = input(1); + x_0_new.D_s_ratio = input(2); + x_0_new.h_f_ratio = input(3); + x_0_new.T_s_ratio = input(4); + x_0_new.F_max = input(5); + x_0_new.D_int = input(6); + x_0_new.w_n = input(7); + [Xs_opt_original, ~, ~, ~, lambda_original, gs] = gradient_optim(x_0_new,p,b,which_objs); + g(:,i) = gs; + Xs_opt_original(8,:) = []; + Xs_opt(:,i) = Xs_opt_original; + lambda.active(:,i) = lambda_original.ineqnonlin; + lambda.lower(:,i) = lambda_original.lower; + lambda.upper(:,i) = lambda_original.upper; + end + residuals2.ineqnonlin = g'; + residuals2.lower = Xs_opt'-b.X_mins'; + residuals2.upper = Xs_opt'-b.X_maxs'; + [~]=constraint_active_plot(residuals2,fval,tol); % save mat file to be read by pareto_bruteforce.m - save('optimization/multiobjective/pareto_search_results',"fval","x","residuals") + save('optimization/multiobjective/pareto_search_results',"fval","x","residuals",... + "lambda") + end function [idx] = constraint_active_plot(residuals,fval,tol) @@ -74,15 +105,58 @@ [~,idx] = sort(fval(:,1)); % order by increasing LCOE figure - subplot 311 + + %subplot 311 + H=subplot(3,1,1, 'Position', [0.195 0.787 0.86 0.16]); spy(lb_active(idx,:)'); - title('Lower Bound Active') + yticks(linspace(1,8,8)); + yticklabels({'WEC Surface Float Outer Diameter', ... + 'Ratio of WEC Surface Float Inner Diameter to Outer Diameter', ... + 'Ratio of WEC Surface Float Height to Outer Diameter', ... + 'Percent of WEC Spar Submergence','Maximum Powertrain Force', ... + 'Powertrain/Controller Damping','Controller Natural Frequency', ... + 'Material'}); + xlabel("Designs, ordered by increasing LCOE", "Fontsize", 9) + grid on - subplot 312 + %subplot 312 + I=subplot(3,1,2, 'Position', [0.195 0.505 0.86 0.16]); spy(ub_active(idx,:)') - title('Upper Bound Active') + yticks(linspace(1,8,8)) + yticklabels({'WEC Surface Float Outer Diameter', ... + 'Ratio of WEC Surface Float Inner Diameter to Outer Diameter', ... + 'Ratio of WEC Surface Float Height to Outer Diameter', ... + 'Percent of WEC Spar Submergence','Maximum Powertrain Force', ... + 'Powertrain/Controller Damping','Controller Natural Frequency', ... + 'Material'}); + xlabel("Designs, ordered by increasing LCOE", "Fontsize", 9) + grid on - subplot 313 + %subplot 313 + J=subplot(3,1,3, 'Position', [0.2 0.08 0.85 0.3]); spy(con_active(idx,:)') - title('Constraint Active') + xlabel("Designs, ordered by increasing LCOE", "Fontsize", 9) + yticks(linspace(1,14,14)); + yticklabels({'Prevent Float Too Heavy', 'Prevent Float Too Light', ... + 'Prevent Spar Too Heavy','Prevent Spar Too Light', ... + 'Metacentric Height','Float Yield','Column Yield','Reaction Plate Yield' ... + 'Column Buckling','Net Generated Power','Minimum Damping Plate Diameter', ... + 'Prevent Float Above Top of the Spar','Prevent LCOE Greater Than Nominal', ... + 'Prevent Irrelevant Max Force'}); + grid on + + improvePlot + + %subplot 311 + set(H, "FontSize",10) + title(H,'Lower Bound Active', "FontSize", 16) + + %subplot 312 + set(I, "FontSize",10) + title(I,'Upper Bound Active', "FontSize", 16) + + %subplot 313 + set(J, "FontSize",9.5) + title(J,'Constraint Active', "FontSize", 16) + end \ No newline at end of file