Anchor.getSize updates to get the design optimization working for SP/DG

- Small fix in project.addMooring related to the "id_part" list
- Took all the hard-coded print statements in capacity_suction (of anchors_famodel) out as there were so many when running Anchor.getSize and calling getAnchorCapacity each iteration
- Anchor.getSize updates to get COBYLA running well for suction piles and D&G piles
- - constraint functions need all the args listed out one by one where the objective function only needs an args parameter...not sure why right now...might be a better way to do this
- - made sure getAnchorCapacity was run in each constraint function (in case some constraints get turned off for different applications)...running getAnchorCapacity is very quick for suction piles (but a little slow for D&G piles)
- - changed the conval for the unity check constraint to make sure the unity check was less than 1 and not greater than 1
- - added a constraint for designing D&G anchors based on the lateral and rotational displacements
- - added a constraint for the design variable bounds, since COBYLA doesn't work well with bounds
- - took out the conFunH and conFunV constraint functions as the unity check constraint appeared to take care of those constraints as well (but this may need to be updated in the future)
- Changed inputs to getFS to the input loads value, rather than using the self.loads value
- Noticed that the return of "getLugForces" wasn't returning what the output was (mudloads)

Author: shousner

famodel/anchors/anchor.py

@@ -623,7 +623,7 @@ def makeEqual_TaTm(mudloads):
         else:
             mudloads['method'] = mudloads['method']
 
-        return(self.loads)
+        return mudloads
 
     def getFS(self, loads=None, acceptance_crit=None):
         '''
@@ -840,7 +840,12 @@ def getSize(self, geom, geomKeys, geomBounds=None, loads=None, minfs={'Ha':1.6,'
         # - - - - Objective and Constraint Functions 
 
         # Define the objective function: Minimize weight of anchor (cost is dependent on weight)
-        def objective(vars,geomKeys, input_loads, fix_zlug):
+        def objective(vars, args):
+
+            geomKeys = args['geomKeys']
+            input_loads = args['input_loads']
+            fix_zlug = args['fix_zlug']
+
             newGeom = dict(zip(geomKeys,vars))
             self.dd['design'].update(newGeom)
             if 'suction' in self.dd['type'] and not fix_zlug:
@@ -851,8 +856,12 @@ def objective(vars,geomKeys, input_loads, fix_zlug):
             return(results['Weight'])
 
         # constraint for suction bucket sizing only. May add more constraints for other anchors in the future...
-        def conFun_LD(vars,LD_con,geomKeys):
-            newGeom = dict(zip(geomKeys,vars))
+        def conFun_LD(vars, geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs):
+            newGeom = dict(zip(geomKeys, vars))
+            self.dd['design'].update(newGeom)
+
+            results = self.getAnchorCapacity(loads=input_loads, plot=False)
+            
             convalA = newGeom['L']/newGeom['D'] - LD_con[0]
             convalB = LD_con[1] - newGeom['L']/newGeom['D']
             conval = min([convalA,convalB])
@@ -863,13 +872,22 @@ def conFun_LD(vars,LD_con,geomKeys):
 
             return(conval)
         # constraint to ensure unity check > 1 for suction buckets
-        def conFun_Suction(vars,input_loads):
-            results = self.getAnchorCapacity(loads=input_loads,plot=False)
-            conval = results['UC'] - 1
+        def conFun_Suction(vars, geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs):
+            results = self.getAnchorCapacity(loads=input_loads, plot=False)
+            #conval = results['UC'] - 1
+            conval = 1 - results['UC']
             # convalB = 1 - results['UC'] 
             return(conval)
+        
+        def conFun_DandG(vars, geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs):
+
+            newGeom = dict(zip(geomKeys, vars))
+            self.dd['design'].update(newGeom)
+            results = self.getAnchorCapacity(loads=input_loads, plot=False)
+
+            return np.array([0.05*newGeom['D'] - results['Lateral displacement'] , 0.25 - results['Rotational displacement']])
 
-        def conFunH(vars,input_loads):
+        def conFunH(vars, geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs):
             # if 'suction' in self.dd['type']:
             #     results = self.getAnchorCapacity(plot=False)
             #     conval = results['UC'] - 1
@@ -885,8 +903,8 @@ def conFunH(vars,input_loads):
                 #             conval = -1*(1-val/minfs[key])
             return(conval)
 
-        def conFunV(vars,minfs,input_loads):
-            FS = self.getFS(input_loads)
+        def conFunV(vars, geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs):
+            FS = self.getFS(loads=input_loads)
             # special case for DEAs
             if minfs['Va'] == 0:
                  conval = 1
@@ -896,6 +914,20 @@ def conFunV(vars,minfs,input_loads):
             # print('FS_V',FS['Va'])
             return(conval)
 
+        def conBounds(vars, geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs):
+
+            newGeom = dict(zip(geomKeys, vars))
+            self.dd['design'].update(newGeom)
+
+            results = self.getAnchorCapacity(loads=input_loads, plot=False)
+
+            bound_L_lower = newGeom['L'] - geomBounds[0][0]
+            bound_L_upper = geomBounds[0][1] - newGeom['L']
+            bound_D_lower = newGeom['D'] - geomBounds[1][0]
+            bound_D_upper = geomBounds[1][1] - newGeom['D']
+
+            return np.array([bound_L_lower, bound_L_upper, bound_D_lower, bound_D_upper])
+        
         # - - - - - Setup & Optimization
         from scipy.optimize import minimize
         from copy import deepcopy
@@ -940,14 +972,20 @@ def conFunV(vars,minfs,input_loads):
         if 'suction' in anchType:
             # bounds = [(1, 7), (5, 50),()]   # Bounds for D and L
             # constraints
-            constraints = [{'type':'ineq','fun':conFun_LD,'args':(LD_con,geomKeys)},
-                           {'type':'ineq','fun':conFun_Suction,'args':(input_loads,)},
-                           {'type':'ineq','fun':conFunH,'args':(input_loads,)},
-                           {'type':'ineq','fun':conFunV,'args':(minfs,input_loads)}]
+            
+            constraints = [{'type':'ineq','fun':conFun_LD,'args':(geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs)},
+                           {'type':'ineq','fun':conFun_Suction,'args':(geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs)},
+                           #{'type':'ineq','fun':conFunH,'args':(geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs)},
+                           #{'type':'ineq','fun':conFunV,'args':(geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs)},
+                           {'type':'ineq','fun':conBounds,'args':(geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs)}]
+        
         else:
             # bounds = None
-            constraints = ({'type':'ineq','fun':conFunH,'args':(input_loads,)},
-                           {'type':'ineq','fun':conFunV,'args':(minfs,input_loads)})
+            constraints = [{'type':'ineq','fun':conFun_LD,'args':(geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs)},
+                           {'type':'ineq','fun':conFun_DandG,'args':(geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs)},
+                           #{'type':'ineq','fun':conFunH,'args':(geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs)},
+                           #{'type':'ineq','fun':conFunV,'args':(geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs)},
+                           {'type':'ineq','fun':conBounds,'args':(geomKeys, input_loads, fix_zlug, LD_con, geomBounds, minfs)}]
 
         # Run the optimization to find sizing that satisfy UC close to 1
         print('optimizing anchor size')
@@ -956,11 +994,11 @@ def conFunV(vars,minfs,input_loads):
             solution = minimize(objective, initial_guess, args=(geomKeys, input_loads, fix_zlug), method="COBYLA",
                                 constraints=constraints, options={'rhobeg':0.1, 'catol':0.001})
         else:
-            solution = minimize(objective, initial_guess, args=(geomKeys, input_loads, fix_zlug), method="COBYLA",
-                                constraints=constraints, bounds=geomBounds, 
+            solution = minimize(objective, initial_guess, args=dict(geomKeys=geomKeys, input_loads=input_loads, fix_zlug=fix_zlug, LD_con=LD_con, geomBounds=geomBounds, minfs=minfs), method="COBYLA",
+                                constraints=constraints,
                                 options={'rhobeg':0.1, 'catol':0.001})
 
-        FS, acceptance, FSdiff = self.getFS(acceptance_crit=minfs)
+        FS, acceptance, FSdiff = self.getFS(loads=input_loads, acceptance_crit=minfs)
 
         # adjust starting value if you're far off from the acceptance criteria (in either direction)
         if FSdiff_max:
@@ -994,7 +1032,7 @@ def conFunV(vars,minfs,input_loads):
                                         constraints=constraints, bounds=geomBounds, 
                                         options={'rhobeg':0.1, 'catol':0.001,'maxiter':400})
                 # re-determine FS and diff from minFS
-                FS, acceptance, FSdiff = self.getFS(acceptance_crit=minfs)  
+                FS, acceptance, FSdiff = self.getFS(loads=input_loads, acceptance_crit=minfs)  
                 count += 1
 
         # Extract the optimized values of geometry

famodel/anchors/anchors_famodel/capacity_suction.py

@@ -83,7 +83,7 @@ def calc_delta(Dr_val):
     if soil_type == 'clay':
         # Definitions for cohesive soils
         Nc = min (6.2*(1 + 0.34*np.arctan(lambdap)),9)   # End-bearing capacity factor
-        ez = (Su0*L**2/2 + k*L**3/3)/(Su0*L + k*L**2/2); print(ez)
+        ez = (Su0*L**2/2 + k*L**3/3)/(Su0*L + k*L**2/2)#; print(ez)
         Np_fixed = 10.25; Np_free = 4                    # From Np vs L/D chart from CAISSON_VHM
         Su_av_L = Su0 + k*zlug                           # Undrained shear strength values (average) 
         Su_tip = Su0 + k*L                               # Undrained shear strength values (tip)
@@ -126,15 +126,15 @@ def f(Hmax):
         # "Leaking"        
         Vmax3 = (PileWeight(L, D, t, rhows) + PileSurface(L, D)*alphastar*Su_av_L + SoilWeight(L, D, t, gamma))                  
         Vmax = min(Vmax1, Vmax2, Vmax3)
-        
+        '''
         print("\n--- Parameter-Based Version ---")
         print(f"Su_av_L        = {Su_av_L:.3f} kPa")
         print(f"sigma'_v(zlug) = {sigma_v_eff:.3f} kPa")
         print(f"psi_val        = {psi_val:.3f}")
         print(f"alpha (API)    = {alpha:.3f}")
         print(f"Hmax           = {Hmax[0]:.2f} kN")
         print(f"Vmax           = {Vmax:.2f} kN")
-
+        '''
 
     elif soil_type == 'sand':
         # Definition for non-cohesive soils
@@ -168,22 +168,22 @@ def f(Hmax):
             # return np.e**(-depth) - 1 + depth
         # Ze = D/(4*7); Zi = D/(4*5)    
         # Vmax = 7*gamma*Ze**2*y(L/Ze)*PileSurface(L, D)/L + 5*gamma*Zi**2*y(L/Zi)*PileSurface(L,(D - 2*t))/L
-        
+        '''
         print("\n--- Parameter-Based (Sand) ---")
         print(f"phi        = {phi:.2f} deg")
         print(f"gamma      = {gamma:.2f} kN/m3")
         print(f"deltastar  = {deltastar:.2f} -")
         print(f"sigma_av_L = {sigma_av_L:.2f} kN")
         print(f"sigma_tip  = {sigma_tip:.2f} kN")
-   
+        '''
     # Pile weight (inc. stiffening plus vent) assessed as a factor
     Wp = 1.10*PileWeight(L, D, t, (rhows + rhow)) 
     # Submerged weight of the soil plug
     Ws = SoilWeight(L, D, t, gamma)
 
     # Capacity envelope
     aVH = 0.5 + lambdap; bVH = 4.5 + lambdap/3 
-    print('Env. exp = ' +str(aVH)+'   '+str(bVH))
+    #print('Env. exp = ' +str(aVH)+'   '+str(bVH))
     UC = (H/Hmax)**aVH + (V/Vmax)**bVH      
     x = np.cos(np.linspace (0, np.pi/2, 100))
     y = (1 - x**bVH)**(1/aVH)

famodel/project.py

@@ -1564,8 +1564,7 @@ def addMooring(self, id=None, endA=None, endB=None, heading=0, dd={},
             elif len(id_part)==1:
                 # anchored line
                 n_existing_moor = len(self.platformList[id_part[0]].getMoorings())
-
-                id = str(id_part)+alph[n_existing_moor]
+                id = str(id_part[0])+alph[n_existing_moor]
             else:
                 id = 'moor'+str(len(self.mooringList))