Local param sensitivity

mdocean/optimization/find_nominal_inputs.m

@@ -15,7 +15,7 @@
         X = [b.X_noms; 1];
         X(5) = F_max_nom;
         [LCOE, P_var, ~, g] = simulation(X,p)
-        [feasible,failed] = is_feasible(g, X, p, b)
+        [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)'})

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, lambdas, grads, hesses] = gradient_optim(x0_input,p,b,which_objs)
 
 if nargin == 0
     % set default parameters if function is run without input
@@ -43,21 +43,22 @@
                                 'PlotFcn',plotfn,...
                                 'MaxIterations',8,...
                                 'FunValCheck','on',...
-                                'ConstraintTolerance',1e-5);%,...
+                                'ConstraintTolerance',1e-5,...
+                                'FiniteDifferenceStepSize',1e-4);%,...
                                 %'SpecifyObjectiveGradient',true,...
                                 %'SpecifyConstraintGradient',true); % would require ALL constraints to be AD-supported
 
 % 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);
+    [Xs_opt, objs_opt, flags, probs, lambdas, grads, hesses] = 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, lambdas, grads, hesses] = optimize_both_objectives(X,p,b,x0_input,opts,ploton,which_objs)
 
     num_constraints = length(b.constraint_names);
     num_objectives = length(which_objs);
@@ -70,34 +71,26 @@
     obj_names = {'LCOE','P_var'};
     probs = cell([1 2]); 
 
+    % allocate outputs
     Xs_opt = zeros(length(X),num_objectives);
     objs_opt = zeros(1,num_objectives);
     flags = zeros(1,num_objectives);
+    grads = zeros(length(X)-1,num_objectives);
+    hesses = zeros(length(X)-1,length(X)-1,num_objectives);
 
-    % add constraints
+    % add nonlinear constraints
     prob = optimproblem();
     for i = 1:num_constraints
         name = b.constraint_names{i};
         prob.Constraints.(name) = g(i) >= 0;
     end
 
-    D_s   = X(1);     % inner diameter of float (m)
-    D_f   = X(2);     % normalized diameter of spar column (-)
-    T_f_2 = X(3);     % normalized draft of float (-)
-    h_s   = X(4);     % normalized draft of spar (-)
-
-    T_s = D_s * p.T_s_over_D_s;
-    h_f = T_f_2 / p.T_f_2_over_h_f;
-    D_d = p.D_d_over_D_s * D_s;
-    D_f_in = D_s * p.D_f_in_over_D_s;
-
-    MEEM = pi*p.harmonics / (p.besseli_argmax*2);
-    prob.Constraints.linear_spar_natural_freq = D_d >= p.D_d_min;
-    prob.Constraints.linear_float_spar_diam = D_f_in <= D_f - .01;
-    prob.Constraints.linear_float_spar_draft = T_f_2 <= T_s - .01;
-    prob.Constraints.linear_float_spar_tops = h_s - T_s >= h_f - T_f_2 + .01;
-    prob.Constraints.linear_float_seafloor = p.h - T_f_2 >= MEEM * D_f; % M
-    prob.Constraints.linear_spar_seafloor = p.h - T_s >= MEEM * D_s; % N
+    % add linear constraints
+    [A_lin,b_lin] = lin_ineq_constraints(p);
+    for i=1:length(b_lin)
+        name = b.lin_constraint_names{i};
+        prob.Constraints.(name) = A_lin(i,:)*X(1:4).' <= b_lin(i);
+    end
 
     % iterate through the two objectives: LCOE and P_var
     for i = 1:num_objectives
@@ -139,6 +132,13 @@
         Xs_opt(:,i) = X_opt;
         objs_opt(i) = obj_opt;
         flags(i) = flag;
+        if i==1
+            lambdas = lambda;
+        else
+            lambdas(i) = lambda;
+        end
+        grads(:,i) = grad;
+        hesses(:,:,i) = hess;
 
         % Post process
         if ploton

mdocean/optimization/run_solver.m

@@ -44,10 +44,18 @@
             % Run scaled optimization problem
             [X_opt,obj_opt,flag,output,lambda,grad,hess] = fmincon(problem_s);
             X_opt = scale .* X_opt;
+            lambda.lower = inv_scale .* lambda.lower;
+            lambda.upper = inv_scale .* lambda.upper;
+            grad = inv_scale .* grad;
+            hess = (inv_scale * inv_scale.') .* hess;
         end
 
         % Rearrange outputs
         X_opt = X_opt(idxs); % reorder elements based on order in autogenerated objective files
+        lambda.lower = lambda.lower(idxs);
+        lambda.upper = lambda.upper(idxs);
+        grad = grad(idxs);
+        hess = hess(idxs,idxs);
 
     else
         [opt_x, obj_opt, flag,output,lambda] = solve(prob,x0,'Options',opts);