symbiotic-engineering · rebeccamccabe · Mar 16, 2024 · Mar 16, 2024 · Mar 16, 2024 · Mar 16, 2024
diff --git a/mdocean/inputs/parameters.m b/mdocean/inputs/parameters.m
@@ -50,6 +50,7 @@
             'LCOE_max', .5,...                  % maximum LCOE ($/kWh)
             'power_max', Inf,...                % maximum power (W)
             'eff_pto',   0.80,...               % PTO efficiency (-)
-            'eff_array', 0.95*0.98 );           % array availability and transmission efficiency (-)
+            'eff_array', 0.95*0.98, ...         % array availability and transmission efficiency (-)
+            'freq_based_optimal_ctrl', false);  % controls based on frequency
 
 end
diff --git a/mdocean/inputs/var_bounds.m b/mdocean/inputs/var_bounds.m
@@ -53,6 +53,9 @@
 b.constraint_names = {'float_too_heavy','float_too_light','spar_too_heavy','spar_too_light',...
     'stability','FOS_float_yield','FOS_col_yield','FOS_plate','FOS_col_buckling',...
     'pos_power','spar_damping','spar_height','LCOE_max','irrelevant_max_force'};
+for i = 15:42
+   b.constraint_names{i} = strcat('prevent_slamming',num2str(i-14));
+end
 
 % modify nominal control inputs to so power and force matches actual
 [F_max_nom, B_p_nom, w_n_nom] = find_nominal_inputs(b, false);

diff --git 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, lambda, gs] = gradient_optim(x0_input,p,b,which_objs)
 
 if nargin == 0
     % set default parameters if function is run without input
@@ -40,17 +40,17 @@
 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]},...
+                                            'OutputSize',{[1,1],[1,1],size(p.JPD),[1, length(b.constraint_names)]},...
                                             'ReuseEvaluation',true,'Analysis','off');%simulation(X, p);
 
     objs = [LCOE P_var];
@@ -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
diff --git 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,133 @@ 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)
+    % I'll change these later
+    for i=1:m
+        for j = 1:n
+            if rem(j,14)==1
+                color = '.r';
+            elseif rem(j,14)==2
+                color = 'or';
+            elseif rem(j,14)==3
+                color = '+r';
+            elseif rem(j,14)==4
+                color = '*r';
+            elseif rem(j,14)==5
+                color = 'xr';
+            elseif rem(j,14)==6
+                color = '-r';
+            elseif rem(j,14)==7
+                color = '|r';
+            elseif rem(j,14)==8
+                color = '+b';
+            elseif rem(j,14)==9
+                color = 'ob';
+            elseif rem(j,14)==10
+                color = '+b';
+            elseif rem(j,14)==11
+                color = '*b';
+            elseif rem(j,14)==12
+                color = 'xb';
+            elseif rem(j,14)==13
+                color = '-b';
+            else
+                color = '|b';
+            end
+            plot(i*100/60,lambdaActive(i,j),color,'MarkerSize',12);
+        end
+    end
+    %plot(pct_angle,lambda_active_sorted)
+    title("Lagrange Multipliers")
+    xlabel("Percent Along the Pareto Front")
+    ylabel("Lagrange Multiplier")
+    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
+    %improvePlot
+
+    % lower bound plot
+    figure
+    hold on
+    [m,n]=size(lambdaLower);
+    % I'll change these later
+    for i=1:m
+        for j = 1:n
+            if rem(j,7)==1
+                color = '.r';
+            elseif rem(j,7)==2
+                color = 'or';
+            elseif rem(j,7)==3
+                color = '+r';
+            elseif rem(j,7)==4
+                color = '*r';
+            elseif rem(j,7)==5
+                color = 'xr';
+            elseif rem(j,7)==6
+                color = '-r';
+            else
+                color = '|r';
+            end
+            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
+            if rem(j,7)==1
+                color = '.r';
+            elseif rem(j,7)==2
+                color = 'or';
+            elseif rem(j,7)==3
+                color = '+r';
+            elseif rem(j,7)==4
+                color = '*r';
+            elseif rem(j,7)==5
+                color = 'xr';
+            elseif rem(j,7)==6
+                color = '-r';
+            else
+                color = '|r';
+            end
+            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')
+        hold off
 end
diff --git a/mdocean/optimization/multiobjective/pareto_search.m b/mdocean/optimization/multiobjective/pareto_search.m
@@ -60,29 +60,103 @@
     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)
     lb_active = abs(residuals.lower) < tol;
     ub_active = abs(residuals.upper) < tol;
-    con_active = abs(residuals.ineqnonlin) < tol;
+    con_active = abs(residuals.ineqnonlin) < tol
 
     [~,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