Added more precise calculation of Cp-optimizing pitch and updated the…

… FBP and VBP tuning procedure to compute everything based on min_pitch rather than assuming fine pitch to be zero

rosco/toolbox/controller.py

@@ -83,7 +83,6 @@ def __init__(self, controller_params):
         self.interp_type = controller_params['interp_type']
 
         # Optional parameters with defaults
-        self.min_pitch          = controller_params['min_pitch']
         self.max_pitch          = controller_params['max_pitch']
         self.vs_minspd          = controller_params['vs_minspd']
         self.ss_vsgain          = controller_params['ss_vsgain']
@@ -106,7 +105,11 @@ def __init__(self, controller_params):
         self.fbp_U = controller_params['VS_FBP_U'] # DBS: Should we set this default based on rated speed?
         self.fbp_P = controller_params['VS_FBP_P']
 
-        #  Optional parameters without defaults
+        # Optional parameters without defaults
+        if 'min_pitch' in controller_params:
+            # Set below if no user input
+            self.min_pitch = controller_params['min_pitch']
+
         if self.VS_FBP > 0:
 
             # Fail if generator torque enabled in Region 3 but pitch control not disabled (may enable these modes to operate together in the future)
@@ -198,7 +201,7 @@ def __init__(self, controller_params):
                 not len(self.U_pc) == len(self.omega_pc) == len(self.zeta_pc):
             raise Exception(
                 'U_pc, omega_pc, and zeta_pc are all list-like and are not of equal length')
-        
+
 
     def tune_controller(self, turbine):
         """
@@ -220,6 +223,9 @@ def tune_controller(self, turbine):
 
         # ------------- Saturation Limits --------------- #
         turbine.max_torque = turbine.rated_torque * self.controller_params['max_torque_factor']
+        if not hasattr(self, 'min_pitch') or not isinstance(self.min_pitch, float): # Set min pitch to optimal if no user input
+            self.min_pitch = turbine.Cp.pitch_opt
+        turbine.min_pitch = self.min_pitch
 
         # -------------Define Operation Points ------------- #
         TSR_rated = rated_rotor_speed*R/turbine.v_rated  # TSR at rated
@@ -243,17 +249,18 @@ def tune_controller(self, turbine):
             if self.fbp_power_mode == 0:
                 P_user_defined *= turbine.rated_power
             # Maximum potential power from MPPT (extending Region 2 power curve to cut-out)
-            P_max = 0.5 * turbine.rho * np.pi*turbine.rotor_radius**2 * turbine.Cp.max * v**3 \
+            Cp_operational = turbine.Cp.interp_surface(self.min_pitch, turbine.TSR_operational) # Compute TSR at controller pitch saturation and operational TSR (may not be optimal)
+            P_max = 0.5 * turbine.rho * np.pi*turbine.rotor_radius**2 * Cp_operational * v**3 \
                 * turbine.GBoxEff/100 * turbine.GenEff/100 # Includes generator efficiency reduction from available inflow power
-            # Take minimum
+            # Take minimum between user input and max possible
             P_op = np.min([P_user_defined, P_max], axis=0)
             # Operation along Cp surface (with fixed pitch)
-            Cp_op = (P_op / P_max) * turbine.Cp.max
+            Cp_op = (P_op / P_max) * Cp_operational
             Cp_op_br = Cp_op[:len(v_below_rated)]
             Cp_op_ar = Cp_op[len(v_below_rated):]
 
             # Identify TSR matching the Cp values (similar to variable pitch angle interpolation below)
-            Cp_FBP = np.ndarray.flatten(turbine.Cp.interp_surface(0, turbine.TSR_initial))     # all Cp values for fine blade pitch (assumed 0)
+            Cp_FBP = np.ndarray.flatten(turbine.Cp.interp_surface(self.min_pitch, turbine.TSR_initial))     # all Cp values for fine blade pitch
             Cp_maxidx = Cp_FBP.argmax()
             # When we depart from Cp_max, our TSR has to fall to either above or below TSR_opt, leading to overspeed and underspeed configurations
             if self.fbp_speed_mode: # Overspeed
@@ -265,21 +272,23 @@ def tune_controller(self, turbine):
                 Cp_op = np.clip(Cp_op, np.min(Cp_FBP[:Cp_maxidx+1]), np.max(Cp_FBP[:Cp_maxidx+1]))            # saturate Cp values to be on Cp surface                                                             # Find maximum Cp value for this TSR
                 f_cp_TSR = interpolate.interp1d(Cp_FBP[:Cp_maxidx+1], turbine.TSR_initial[:Cp_maxidx+1])             # interpolate function for Cp(tsr) values
             TSR_op = f_cp_TSR(Cp_op)
-            TSR_below_rated = TSR_op[:len(v_below_rated)]
+            # Defer to operational TSR for below rated, even if other optimum (should keep Cp the same, but may lead to discontinuities in operating schedule if min_pitch is not optimal)
+            TSR_op[v < turbine.v_rated] = turbine.TSR_operational
+            TSR_below_rated = TSR_op[:len(v_below_rated)] # Should be constant at TSR_operational
             TSR_above_rated = TSR_op[len(v_below_rated):]
 
         # elif self.PC_ControlMode > 0: # If using pitch control in Region 3
         else: # Default here even if pitch control disabled to maintain backwards compatibility
 
-            # separate TSRs by operations regions
-            TSR_below_rated = [min(turbine.TSR_operational, rated_rotor_speed*R/v) for v in v_below_rated] # below rated
+            # TSR setpoints
+            TSR_below_rated = [min(turbine.TSR_operational, rated_rotor_speed*R/v) for v in v_below_rated] # below rated, saturated to not exceed rated rotor speed
             TSR_above_rated = rated_rotor_speed*R/v_above_rated                     # above rated
             TSR_op = np.concatenate((TSR_below_rated, TSR_above_rated))             # operational TSRs
 
             # Find expected operational Cp values
-            Cp_above_rated = turbine.Cp.interp_surface(0,TSR_above_rated[0])     # Cp during rated operation (not optimal). Assumes cut-in bld pitch to be 0
-            Cp_op_br = np.ones(len(v_below_rated)) * turbine.Cp.max              # below rated
+            Cp_above_rated = turbine.Cp.interp_surface(self.min_pitch,TSR_above_rated[0])     # Cp during rated operation (not optimal)
             Cp_op_ar = Cp_above_rated * (TSR_above_rated/TSR_rated)**3           # above rated
+            Cp_op_br = [turbine.Cp.interp_surface(self.min_pitch, TSR) for TSR in TSR_below_rated] # Below rated, reinterpolated for accurate operational Cp based on controller pitch
             Cp_op = np.concatenate((Cp_op_br, Cp_op_ar))                         # operational CPs to linearize around
 
 
@@ -299,11 +308,8 @@ def tune_controller(self, turbine):
 
             if self.VS_FBP > 0: # Fixed blade pitch control in Region 3
 
-                if isinstance(self.min_pitch, float):
-                    pitch_op[i] = self.min_pitch
-                else:
-                    # pitch_op[i] = 0 # Assume zero pitch
-                    pitch_op[i] = turbine.pitch_initial_rad[turbine.Cp.max_ind[1]] # Take optimal pitch from Cp surface
+                # Constant pitch, either user input or optimal pitch from Cp surface
+                pitch_op[i] = self.min_pitch
 
             else: # Variable pitch control in Region 3 (default)
 
@@ -314,12 +320,10 @@ def tune_controller(self, turbine):
                 f_cp_pitch = interpolate.interp1d(Cp_TSR[Cp_maxidx:],pitch_initial_rad[Cp_maxidx:])             # interpolate function for Cp(tsr) values
 
                 # expected operational blade pitch values. Saturates by min_pitch if it exists
-                if v[i] <= turbine.v_rated and isinstance(self.min_pitch, float): # Below rated & defined min_pitch
+                if v[i] <= turbine.v_rated: # Below rated, keep lowest pitch
                     pitch_op[i] = min(self.min_pitch, f_cp_pitch(Cp_op[i]))
-                elif isinstance(self.min_pitch, float):                           # above rated & defined min_pitch
+                else:                       # above rated, keep highest pitch
                     pitch_op[i] = max(self.min_pitch, f_cp_pitch(Cp_op[i]))
-                else:                                                             # no defined minimum pitch schedule
-                    pitch_op[i] = f_cp_pitch(Cp_op[i])     
 
             # Calculate Cp Surface gradients
             dCp_beta[i], dCp_TSR[i] = turbine.Cp.interp_gradient(pitch_op[i],TSR_op[i]) 
@@ -331,10 +335,6 @@ def tune_controller(self, turbine):
             Ct_op[i]    = f_ct(pitch_op[i])
             Ct_op[i]    = np.clip(Ct_op[i], np.min(Ct_TSR), np.max(Ct_TSR))        # saturate Ct values to be on Ct surface
 
-        # Define minimum pitch saturation to be at Cp-maximizing pitch angle if not specifically defined
-        if not isinstance(self.min_pitch, float):
-            self.min_pitch = pitch_op[0]
-
         # Compute generator speed and torque operating schedule
         P_op = 0.5 * turbine.rho * np.pi*turbine.rotor_radius**2 * Cp_op * v**3 * turbine.GBoxEff/100 * turbine.GenEff/100
         if self.VS_FBP == 0: # Saturate between min speed and rated if variable pitch in Region 3
@@ -411,8 +411,10 @@ def tune_controller(self, turbine):
         self.pc_gain_schedule.second_order_PI(self.zeta_pc_U, self.omega_pc_U,A_pc,B_beta[-len(v_above_rated)+1:],linearize=True,v=v_above_rated[1:])        
         self.vs_gain_schedule = ControllerTypes()
         if self.VS_FBP == 0:
+            # Using variable pitch control in Region 3, so generate torque gain schedule only for Region 2
             self.vs_gain_schedule.second_order_PI(self.zeta_vs, self.omega_vs,A_vs,B_tau[0:len(v_below_rated)],linearize=False,v=v_below_rated)
         else:
+            # Using fixed pitch torque control in Region 3, so generate torque gain schedule for Regions 2 and 3
             self.vs_gain_schedule.second_order_PI(self.zeta_vs, self.omega_vs,A,B_tau,linearize=False,v=v)
 
         # -- Find K for Komega_g^2 --

rosco/toolbox/inputs/toolbox_schema.yaml

@@ -312,7 +312,6 @@ properties:
             min_pitch:           
                 description: Minimum pitch angle [rad], {default = 0 degrees}
                 type: number
-                default: 0
                 unit: rad
             vs_minspd:           
                 description: Minimum rotor speed [rad/s], {default = 0 rad/s}

rosco/toolbox/turbine.py

@@ -57,7 +57,6 @@ class Turbine():
     load_from_ccblade
     load_from_txt
     generate_rotperf_fast
-    write_rotor_performance
 
     Parameters:
     -----------
@@ -272,6 +271,8 @@ def load_from_fast(
         # Define operational TSR
         if not self.TSR_operational:
             self.TSR_operational = self.Cp.TSR_opt
+        # Compute operational Cp (may not be optimal if TSR_operational set by user)
+        self.Cp_operational = self.Cp.interp_surface(self.Cp.pitch_opt, self.TSR_operational)
 
         # Pull out some floating-related data
         try:
@@ -291,7 +292,6 @@ def load_from_ccblade(self):
                   Dictionary containing fast model details - defined using from InputReader_OpenFAST (distributed as a part of AeroelasticSE)
 
         '''
-        from wisdem.ccblade.ccblade import CCAirfoil, CCBlade
 
         print('Loading rotor performance data from CC-Blade.')
 
@@ -334,7 +334,6 @@ def load_from_ccblade(self):
         self.Ct_table = Ct 
         self.Cq_table = Cq
 
-
     def generate_rotperf_fast(self, openfast_path, FAST_runDirectory=None, run_BeamDyn=False,
                               debug_level=1, run_type='multi'):
         '''
@@ -622,7 +621,7 @@ class RotorPerformance():
     TSR_initial : array_like (rad)
                     An [n x 1] or [1 x n] array containing tip-speed ratios corresponding to  performance_table 
     '''
-    def __init__(self,performance_table, pitch_initial_rad, TSR_initial):
+    def __init__(self, performance_table, pitch_initial_rad, TSR_initial):
 
         # Store performance data tables
         self.performance_table = performance_table          # Table containing rotor performance data, i.e. Cp, Ct, Cq
@@ -635,27 +634,39 @@ def __init__(self,performance_table, pitch_initial_rad, TSR_initial):
         # "Optimal" below rated TSR and blade pitch (for Cp) - note this may be limited by resolution of Cp-surface
         self.max = np.amax(performance_table)
         self.max_ind = np.where(performance_table == np.amax(performance_table))
-        self.pitch_opt = pitch_initial_rad[self.max_ind[1]]
-        # --- Find TSR ---
-        # Make finer mesh for Tip speed ratios at "optimal" blade pitch angle, do a simple lookup. 
-        #       -- nja: this seems to work a little better than interpolating
 
-        # Find the 1D performance table when pitch is at the maximum part of the Cx surface:
-        performance_beta_max = np.ndarray.flatten(performance_table[:,self.max_ind[1][-1]]) # performance metric at the last maximizing pitch angle
-
-        # If there is more than one max pitch angle:
+        # Throw a warning if there is more than one max pitch angle
         if len(self.max_ind[1]) > 1:
             print('rosco.toolbox Warning: repeated maximum values in a performance table and the last one @ pitch = {} rad. was taken...'.format(self.pitch_opt[-1]))
 
-        # Find TSR that maximizes Cx at fine pitch
-        # - TSR to satisfy: max( Cx(TSR, \beta_fine) ) = TSR_opt
-        TSR_fine_ind = np.linspace(TSR_initial[0],TSR_initial[-1],int(TSR_initial[-1] - TSR_initial[0])*100) # Range of TSRs to interpolate accross
-        f_TSR = interpolate.interp1d(TSR_initial,TSR_initial,bounds_error='False',kind='quadratic')          # interpolate function for Cp(tsr) values
-        TSR_fine = f_TSR(TSR_fine_ind) # TSRs at fine pitch
-        f_performance = interpolate.interp1d(TSR_initial,performance_beta_max,bounds_error='False',kind='quadratic')    # interpolate function for Cx(tsr) values
-        performance_fine = f_performance(TSR_fine_ind) # Cx values at fine pitch
-        performance_max_ind = np.where(performance_fine == np.max(performance_fine)) # Find max performance at fine pitch
-        self.TSR_opt = float(TSR_fine[performance_max_ind[0]][0])  # TSR to maximize Cx at fine pitch
+        # Refine optimal point using smooth interpolation -- gives more precise estimates than coarse grid
+        # Select region around maximizing TSR and blade pitch
+        TSR_max_ind = self.max_ind[0][-1]
+        pitch_max_ind = self.max_ind[1][-1]
+        ind_search_num = 3
+        TSR_search_ind = [np.max([0, TSR_max_ind - ind_search_num]), np.min([TSR_max_ind + ind_search_num, len(self.TSR_initial)-1])]
+        pitch_search_ind = [np.max([0, pitch_max_ind - ind_search_num]), np.min([pitch_max_ind + ind_search_num, len(self.pitch_initial_rad)-1])]
+        TSR_search_range = range(TSR_search_ind[0], TSR_search_ind[1]+1)
+        pitch_search_range = range(pitch_search_ind[0], pitch_search_ind[1]+1)
+
+        print('rosco.toolbox: Refining maximum performance estimate...')
+
+        # Generate finer mesh in the proximity of maximizer
+        fine_mesh_scale = 20
+        TSR_fine = np.linspace(self.TSR_initial[TSR_search_ind[0]], self.TSR_initial[TSR_search_ind[1]], np.diff(TSR_search_ind)[0] * fine_mesh_scale)
+        pitch_fine = np.linspace(self.pitch_initial_rad[pitch_search_ind[0]], self.pitch_initial_rad[pitch_search_ind[1]], np.diff(pitch_search_ind)[0] * fine_mesh_scale)
+        pitch_fine_mesh, TSR_fine_mesh = np.meshgrid(pitch_fine, TSR_fine) # Construct mesh grids of fine data to interpolate performance surface
+        performance_fine = interpolate.interpn([TSR_initial[TSR_search_range], pitch_initial_rad[pitch_search_range]], \
+                                               self.performance_table[TSR_search_range, :][:, pitch_search_range], \
+                                               np.array([TSR_fine_mesh, pitch_fine_mesh]).T, \
+                                               bounds_error='False', method='cubic') # Cubic spline interpolation over finer mesh data
+
+        # Save optimal performance, TSR, and pitch
+        self.performance_opt = np.max(performance_fine) # Maximal Cx on fine grid
+        performance_max_ind = np.where(performance_fine == self.performance_opt) # Find maximizer
+        self.TSR_opt = float(TSR_fine[performance_max_ind[0][0]]) # If multiple maximizers, pick lowest TSR
+        self.pitch_opt = float(pitch_fine[performance_max_ind[1][-1]]) # If multiple maximizers, pick highest pitch
+
 
     def interp_surface(self,pitch,TSR):
         '''