Merge remote-tracking branch 'upstream/develop' into weis

.github/workflows/CI_rosco-toolbox.yml

@@ -96,7 +96,7 @@ jobs:
       # Install OpenFAST
       - name: Install OpenFAST
         run: |
-          conda install openfast
+          conda install openfast==2.5.0
           
       # Run examples
       - name: Run examples

README.md

@@ -39,17 +39,15 @@ For LaTeX users:
 ```
 If the ROSCO generic tuning theory and implementation played a roll in your research, please cite the following paper
 ```
-@inproceedings{Abbas_WindTech2019,
-	doi = {10.1088/1742-6596/1452/1/012002},
-	url = {https://doi.org/10.1088%2F1742-6596%2F1452%2F1%2F012002},
-	year = 2020,
-	month = {jan},
-	publisher = {{IOP} Publishing},
-	volume = {1452},
-	pages = {012002},
-	author = {Nikhar J. Abbas and Alan Wright and Lucy Pao},
-	title = {An Update to the National Renewable Energy Laboratory Baseline Wind Turbine Controller},
-	journal = {Journal of Physics: Conference Series}
+@Article{wes-2021-19,
+AUTHOR = {Abbas, N. and Zalkind, D. and Pao, L. and Wright, A.},
+TITLE = {A Reference Open-Source Controller for Fixed and Floating Offshore Wind Turbines},
+JOURNAL = {Wind Energy Science Discussions},
+VOLUME = {2021},
+YEAR = {2021},
+PAGES = {1--33},
+URL = {https://wes.copernicus.org/preprints/wes-2021-19/},
+DOI = {10.5194/wes-2021-19}
 }
 ```
 Additionally, if you have extensively used the [ROSCO](https://github.com/NREL/ROSCO) controller or [WISDEM](https://github.com/wisdem/wisdem), please cite them accordingly. 

ROSCO_toolbox/control_interface.py

@@ -9,7 +9,7 @@
 # CONDITIONS OF ANY KIND, either express or implied. See the License for the
 # specific language governing permissions and limitations under the License.
 
-from ctypes import byref, cdll, c_int, POINTER, c_float, c_char_p, c_double, create_string_buffer, c_int32
+from ctypes import byref, cdll, c_int, POINTER, c_float, c_char_p, c_double, create_string_buffer, c_int32, c_void_p
 import numpy as np
 from numpy.ctypeslib import ndpointer
 
@@ -35,19 +35,31 @@ class ControllerInterface():
 
     """
 
-    def __init__(self, lib_name, param_filename='DISCON.IN'):
+    def __init__(self, lib_name, param_filename='DISCON.IN', **kwargs):
         """
         Setup the interface
         """
         self.lib_name = lib_name
         self.param_name = param_filename
 
-        # Temp fixed parameters
+        # Set default parameters
         # PARAMETERS
         self.DT = 0.1
         self.num_blade = 3
         self.char_buffer = 500
         self.avr_size = 500
+        self.sim_name = 'simDEBUG'
+
+        # Set kwargs, like DT
+        for (k, w) in kwargs.items():
+            try:
+                setattr(self, k, w)
+            except:
+                pass
+
+        self.init_discon()
+
+    def init_discon(self):
 
         # Initialize
         self.pitch = 0
@@ -56,10 +68,12 @@ def __init__(self, lib_name, param_filename='DISCON.IN'):
         self.discon = cdll.LoadLibrary(self.lib_name)
         self.avrSWAP = np.zeros(self.avr_size)
 
-        # Define some avrSWAP parameters
+        # Define some avrSWAP parameters, NOTE: avrSWAP indices are offset by -1 from Fortran
         self.avrSWAP[2] = self.DT
         self.avrSWAP[60] = self.num_blade
-        self.avrSWAP[20] = 1 # HARD CODE initial rot speed = 1 rad/s
+        self.avrSWAP[19] = 1.0 # HARD CODE initial gen speed = 1 rad/s
+        self.avrSWAP[20] = 1.0 # HARD CODE initial rot speed = 1 rad/s
+        self.avrSWAP[82] = 0 # HARD CODE initial nacIMU = 0
         self.avrSWAP[26] = 10 # HARD CODE initial wind speed = 10 m/s
 
 
@@ -76,7 +90,7 @@ def __init__(self, lib_name, param_filename='DISCON.IN'):
         self.aviFAIL = c_int32() # 1
         self.accINFILE = self.param_name.encode('utf-8')
         # self.avcOUTNAME = create_string_buffer(1000) # 'DEMO'.encode('utf-8')
-        self.avcOUTNAME = 'simDEBUG.dbg'.encode('utf-8')
+        self.avcOUTNAME = (self.sim_name + '.RO.dbg').encode('utf-8')
         self.avcMSG = create_string_buffer(1000)
         self.discon.DISCON.argtypes = [POINTER(c_float), POINTER(c_int32), c_char_p, c_char_p, c_char_p] # (all defined by ctypes)
 
@@ -103,7 +117,7 @@ def call_discon(self):
         self.avrSWAP = data
 
 
-    def call_controller(self,t,dt,pitch,torque,genspeed,geneff,rotspeed,ws):
+    def call_controller(self,t,dt,pitch,torque,genspeed,geneff,rotspeed,ws,NacIMU_FA_Acc=0):
         '''
         Runs the controller. Passes current turbine state to the controller, and returns control inputs back
         
@@ -123,6 +137,9 @@ def call_controller(self,t,dt,pitch,torque,genspeed,geneff,rotspeed,ws):
                   rotor speed, (rad/s)
         ws: float
             wind speed, (m/s)
+        NacIMU_FA_Acc : float
+            nacelle IMU accel. in the nodding dir. , (m/s**2)
+            default to 0 (fixed-bottom, simple 1-DOF sim does not include it, but OpenFAST linearizations do)
         '''
 
         # Add states to avr
@@ -136,6 +153,7 @@ def call_controller(self,t,dt,pitch,torque,genspeed,geneff,rotspeed,ws):
         self.avrSWAP[19] = genspeed
         self.avrSWAP[20] = rotspeed
         self.avrSWAP[26] = ws
+        self.avrSWAP[82] = NacIMU_FA_Acc
 
         # call controller
         self.call_discon()
@@ -152,3 +170,17 @@ def show_control_values(self):
         '''
         print('Pitch',self.pitch)
         print('Torque',self.torque)
+
+    def kill_discon(self):
+        '''
+        Unload the dylib from memory: https://stackoverflow.com/questions/359498/how-can-i-unload-a-dll-using-ctypes-in-python
+        This is unix-specific, but there seems to be a windows solution as well
+        '''
+
+        print('Shutting down {}'.format(self.lib_name))
+        handle = self.discon._handle
+        dlclose_func = self.discon.dlclose
+        dlclose_func.argtypes = [c_void_p]
+
+        del self.discon
+        dlclose_func(handle)

ROSCO_toolbox/controller.py

@@ -135,6 +135,24 @@ def __init__(self, controller_params):
             else:
                 self.flp_maxpit = 0.0
 
+        # Gain scheduling/indexing
+        # Number of wind speed breakpoints, default = 60
+        if 'WS_GS_n' in controller_params:
+            self.WS_GS_n = controller_params['WS_GS_n']
+        else:
+            self.WS_GS_n = 60
+
+        # Number of pitch control breakpoints, default = 30
+        if 'PC_GS_n' in controller_params:
+            self.PC_GS_n = controller_params['PC_GS_n']
+        else:
+            self.PC_GS_n = 30
+
+        if self.WS_GS_n <= self.PC_GS_n:
+            raise Exception('Number of WS breakpoints is not greater than pitch control breakpoints')
+
+
+
     def tune_controller(self, turbine):
         """
         Given a turbine model, tune a controller based on the NREL generic controller tuning process
@@ -158,8 +176,9 @@ def tune_controller(self, turbine):
         TSR_rated = rated_rotor_speed*R/turbine.v_rated  # TSR at rated
 
         # separate wind speeds by operation regions
-        v_below_rated = np.linspace(turbine.v_min,turbine.v_rated, num=30)[:-1]             # below rated
-        v_above_rated = np.linspace(turbine.v_rated,turbine.v_max, num=30)             # above rated
+        # add one to above rated because we don't use rated in the pitch control gain scheduling
+        v_below_rated = np.linspace(turbine.v_min,turbine.v_rated, num=self.WS_GS_n-self.PC_GS_n)[:-1]             # below rated
+        v_above_rated = np.linspace(turbine.v_rated,turbine.v_max, num=self.PC_GS_n+1)             # above rated
         v = np.concatenate((v_below_rated, v_above_rated))
 
         # separate TSRs by operations regions
@@ -248,7 +267,7 @@ def tune_controller(self, turbine):
 
         # -- Find gain schedule --
         self.pc_gain_schedule = ControllerTypes()
-        self.pc_gain_schedule.second_order_PI(self.zeta_pc, self.omega_pc,A_pc,B_beta[-len(v_above_rated)+1:],linearize=True,v=v_above_rated[1:])
+        self.pc_gain_schedule.second_order_PI(self.zeta_pc, self.omega_pc,A_pc,B_beta[-len(v_above_rated)+1:],linearize=True,v=v_above_rated[1:])        
         self.vs_gain_schedule = ControllerTypes()
         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)
 

ROSCO_toolbox/utilities.py

@@ -131,7 +131,7 @@ def write_DISCON(turbine, controller, param_file='DISCON.IN', txt_filename='Cp_C
     file.write('!------- WIND SPEED ESTIMATOR ---------------------------------------------\n')
     file.write('{:<13.3f}       ! WE_BladeRadius	- Blade length (distance from hub center to blade tip), [m]\n'.format(turbine.rotor_radius))
     file.write('{:<11d}         ! WE_CP_n			- Amount of parameters in the Cp array\n'.format(1))
-    file.write(          '{}    ! WE_CP - Parameters that define the parameterized CP(lambda) function\n'.format(''.join('{:<2.1f} '.format(0.0) for i in range(4))))
+    file.write('{}              ! WE_CP - Parameters that define the parameterized CP(lambda) function\n'.format(''.join('{:<2.1f} '.format(0.0) for i in range(1))))
     file.write('{:<13.1f}		! WE_Gamma			- Adaption gain of the wind speed estimator algorithm [m/rad]\n'.format(0.0))
     file.write('{:<13.1f}       ! WE_GearboxRatio	- Gearbox ratio [>=1],  [-]\n'.format(turbine.Ng))
     file.write('{:<014.5f}      ! WE_Jtot			- Total drivetrain inertia, including blades, hub and casted generator inertia to LSS, [kg m^2]\n'.format(turbine.J))
@@ -157,7 +157,7 @@ def write_DISCON(turbine, controller, param_file='DISCON.IN', txt_filename='Cp_C
     file.write('\n')
     file.write('!------- TOWER FORE-AFT DAMPING -------------------------------------------\n')
     file.write('{:<11d}         ! FA_KI				- Integral gain for the fore-aft tower damper controller, -1 = off / >0 = on [rad s/m] - !NJA - Make this a flag\n'.format(-1))
-    file.write('{:<13.1f}       ! FA_HPF_CornerFreq	- Corner frequency (-3dB point) in the high-pass filter on the fore-aft acceleration signal [rad/s]\n'.format(0.0))
+    file.write('{:<13.1f}       ! FA_HPFCornerFreq	- Corner frequency (-3dB point) in the high-pass filter on the fore-aft acceleration signal [rad/s]\n'.format(0.0))
     file.write('{:<13.1f}       ! FA_IntSat			- Integrator saturation (maximum signal amplitude contribution to pitch from FA damper), [rad]\n'.format(0.0))
     file.write('\n')
     file.write('!------- MINIMUM PITCH SATURATION -------------------------------------------\n')