Add omega vs. windspeed functionality

Examples/example_04.py

@@ -50,14 +50,21 @@
 write_DISCON(turbine,controller,param_file=param_file, txt_filename=os.path.join(tune_dir,path_params['rotor_performance_filename']))
 
 # Plot gain schedule
-fig, ax = plt.subplots(1,2,constrained_layout=True)
-ax[0].plot(controller.v[len(controller.v_below_rated)+1:], controller.pc_gain_schedule.Kp)
-ax[0].set_xlabel('Wind Speed')
-ax[0].set_ylabel('Proportional Gain')
+fig, ax = plt.subplots(2,2,constrained_layout=True,sharex=True)
+ax = ax.flatten()
+ax[0].plot(controller.v[len(controller.v_below_rated)+1:], controller.omega_pc_U)
+ax[0].set_ylabel('omega_pc')
 
-ax[1].plot(controller.v[len(controller.v_below_rated)+1:], controller.pc_gain_schedule.Ki)
-ax[1].set_xlabel('Wind Speed')
-ax[1].set_ylabel('Integral Gain')
+ax[1].plot(controller.v[len(controller.v_below_rated)+1:], controller.zeta_pc_U)
+ax[1].set_ylabel('zeta_pc')
+
+ax[2].plot(controller.v[len(controller.v_below_rated)+1:], controller.pc_gain_schedule.Kp)
+ax[2].set_xlabel('Wind Speed')
+ax[2].set_ylabel('Proportional Gain')
+
+ax[3].plot(controller.v[len(controller.v_below_rated)+1:], controller.pc_gain_schedule.Ki)
+ax[3].set_xlabel('Wind Speed')
+ax[3].set_ylabel('Integral Gain')
 
 plt.suptitle('Pitch Controller Gains')
 

Examples/example_06.py

@@ -13,6 +13,7 @@
 # Python Modules
 import yaml
 import os
+import matplotlib.pyplot as plt
 # ROSCO toolbox modules 
 from ROSCO_toolbox import controller as ROSCO_controller
 from ROSCO_toolbox import turbine as ROSCO_turbine
@@ -47,6 +48,30 @@
 param_file = os.path.join(this_dir,'DISCON.IN')   # This must be named DISCON.IN to be seen by the compiled controller binary. 
 write_DISCON(turbine,controller,param_file=param_file, txt_filename=path_params['rotor_performance_filename'])
 
+# Plot gain schedule
+fig, ax = plt.subplots(2,2,constrained_layout=True,sharex=True)
+ax = ax.flatten()
+ax[0].plot(controller.v[len(controller.v_below_rated)+1:], controller.omega_pc_U)
+ax[0].set_ylabel('omega_pc')
+
+ax[1].plot(controller.v[len(controller.v_below_rated)+1:], controller.zeta_pc_U)
+ax[1].set_ylabel('zeta_pc')
+
+ax[2].plot(controller.v[len(controller.v_below_rated)+1:], controller.pc_gain_schedule.Kp)
+ax[2].set_xlabel('Wind Speed')
+ax[2].set_ylabel('Proportional Gain')
+
+ax[3].plot(controller.v[len(controller.v_below_rated)+1:], controller.pc_gain_schedule.Ki)
+ax[3].set_xlabel('Wind Speed')
+ax[3].set_ylabel('Integral Gain')
+
+plt.suptitle('Pitch Controller Gains')
+
+if False:
+  plt.show()
+else:
+  plt.savefig(os.path.join(example_out_dir,'06_GainSched.png'))
+
 # Run OpenFAST
 # --- May need to change fastcall if you use a non-standard command to call openfast
 fastcall = 'openfast'

ROSCO_toolbox/controller.py

@@ -63,6 +63,7 @@ def __init__(self, controller_params):
         self.Flp_Mode           = controller_params['Flp_Mode']
 
         # Necessary parameters
+        self.U_pc               = controller_params['U_pc']
         self.zeta_pc            = controller_params['zeta_pc']
         self.omega_pc           = controller_params['omega_pc']
         self.zeta_vs            = controller_params['zeta_vs']
@@ -225,9 +226,16 @@ def tune_controller(self, turbine):
         A_pc = A[-len(v_above_rated)+1:]     # above rated
         B_tau = B_tau * np.ones(len(v))
 
+        # Resample omega_ and zeta_pc at above rated wind speeds
+        if len(self.U_pc) == len(self.omega_pc) == len(self.zeta_pc):
+            self.omega_pc_U = multi_sigma(v_above_rated[1:],self.U_pc,self.omega_pc)
+            self.zeta_pc_U  = multi_sigma(v_above_rated[1:],self.U_pc,self.zeta_pc)
+        else:
+            raise Exception('ROSCO_toolbox: The lengths of U_pc, omega_pc, and zeta_pc must be equal')
+
         # -- 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_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()
         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)
 
@@ -517,10 +525,10 @@ def second_order_PI(self,zeta,om_n,A,B,linearize=False,v=None):
 
         Parameters:
         -----------
-        zeta : int (-)
-               Desired damping ratio 
-        om_n : int (rad/s)
-               Desired natural frequency 
+        zeta : list of floats (-)
+               Desired damping ratio with breakpoints at v
+        om_n : list of floats (rad/s)
+               Desired natural frequency with breakpoints at v
         A : array_like (1/s)
             Plant poles (state transition matrix)
         B : array_like (varies)
@@ -540,3 +548,69 @@ def second_order_PI(self,zeta,om_n,A,B,linearize=False,v=None):
         # Calculate gain schedule
         self.Kp = 1/B * (2*zeta*om_n + A)
         self.Ki = om_n**2/B           
+
+
+# helper functions
+
+def sigma(tt,t0,t1,y0=0,y1=1):
+    ''' 
+    generates timeseries for a smooth transition from y0 to y1 from x0 to x1
+
+    inputs: tt - time indices
+            t0 - start time
+            t1 - end time
+            y0 - start output
+            y1 - end output
+
+    outputs: yy - output timeseries corresponding to tt
+    '''
+
+    a3 = 2/(t0-t1)**3
+    a2 = -3*(t0+t1)/(t0-t1)**3
+    a1 = 6*t1*t0/(t0-t1)**3
+    a0 = (t0-3*t1)*t0**2/(t0-t1)**3
+
+    a = np.array([a3,a2,a1,a0])  
+
+    T = np.vander(tt,N=4)       # vandermonde matrix
+
+    ss = T @ a.T                # base sigma
+
+    yy = (y1-y0) * ss + y0      # scale and offset
+
+    return yy
+
+
+def multi_sigma(xx,x_bp,y_bp):
+    '''
+    Make a sigma interpolation with multiple breakpoints
+
+    Parameters:
+    -----------
+    xx : list of floats (-)
+            new sample points
+    x_bp : list of floats (-)
+            breakpoints
+    y_bp : list of floats (-)
+            function value at breakpoints
+
+    '''
+    yy = np.zeros_like(xx)
+
+    # interpolate sigma functions between all breakpoints
+    for i_sigma in range(0,len(x_bp)-1):
+        ind_i       = (xx >= x_bp[i_sigma]) & (xx < x_bp[i_sigma+1])
+        xx_i        = xx[ind_i]
+        yy_i        = sigma(xx_i,x_bp[i_sigma],x_bp[i_sigma+1],y0=y_bp[i_sigma],y1=y_bp[i_sigma+1])
+        yy[ind_i]   = yy_i
+
+    # add first and last values to beginning and end
+    yy[xx<x_bp[0]]      = y_bp[0]
+    yy[xx>x_bp[-1]]     = y_bp[-1]
+
+    if False:  # debug plot
+        import matplotlib.pyplot as plt
+        plt.plot(xx,yy)
+        plt.show()
+
+    return yy

ROSCO_toolbox/inputs/toolbox_schema.yaml

@@ -187,16 +187,28 @@ properties:
                 maximum: 2
                 default: 0
                 description: Active Power Control Mode (0- no active power control 1- constant active power control, 2- open loop power vs time, 3- open loop power vs. wind speed)
+            U_pc:                       
+                type: array
+                description: List of wind speeds to schedule pitch control zeta and omega 
+                unit: m/s
+                items:
+                    type: number
+                    minimum: 0
+                    uniqueItems: True
             zeta_pc:                       
-                type: number
-                minimum: 0
-                description: Pitch controller desired damping ratio [-]
+                type: array
+                description: List of pitch controller desired damping ratio at U_pc [-]
                 unit: none
+                items:
+                    type: number
+                    minimum: 0
             omega_pc:                     
-                type: number
-                minimum: 0
-                description: Pitch controller desired natural frequency [rad/s]
+                type: array
+                description: List of pitch controller desired natural frequency at U_pc [rad/s]
                 unit: rad/s
+                items:
+                    type: number
+                    minimum: 0
             zeta_vs:                        
                 type: number
                 minimum: 0

Tune_Cases/BAR.yaml

@@ -42,8 +42,9 @@ controller_params:
   Fl_Mode:            0                       # Floating specific feedback mode {0: no nacelle velocity feedback, 1: nacelle velocity feedback}
   Flp_Mode:           2                       # Flap control mode {0: no flap control, 1: steady state flap angle, 2: Proportional flap control}
   # Controller parameters 
-  zeta_pc:            1.2                     # Pitch controller desired damping ratio [-]
-  omega_pc:           0.1                     # Pitch controller desired natural frequency [rad/s]
+  U_pc:               [12]
+  zeta_pc:            [1.2]                     # Pitch controller desired damping ratio [-]
+  omega_pc:           [0.1]                     # Pitch controller desired natural frequency [rad/s]
   zeta_vs:            0.5                     # Torque controller desired damping ratio [-]
   omega_vs:           0.1                     # Torque controller desired natural frequency [rad/s]
   # Only needed if Flp_Mode > 0

Tune_Cases/DTU10MW.yaml

@@ -45,8 +45,9 @@ controller_params:
   # twr_freq:           # None                # Tower natural frequency [rad/s]
   # ptfm_freq:          # None                # Platform natural frequency [rad/s]         
 # Controller parameters
-  zeta_pc:            0.8                   # Pitch controller desired damping ratio [-]
-  omega_pc:           0.3                   # Pitch controller desired natural frequency [rad/s]
+  U_pc:               [12]
+  zeta_pc:            [0.8]                  # Pitch controller desired damping ratio [-]
+  omega_pc:           [0.3]                   # Pitch controller desired natural frequency [rad/s]
   zeta_vs:            0.7                   # Torque controller desired damping ratio [-]
   omega_vs:           0.2                   # Torque controller desired natural frequency [rad/s]
   # Optional - these can be defined, but do not need to be

Tune_Cases/IEA15MW.yaml

@@ -39,8 +39,9 @@ controller_params:
   Fl_Mode:            1                     # Floating specific feedback mode {0: no nacelle velocity feedback, 1: nacelle velocity feedback}
   Flp_Mode:           0                     # Flap control mode {0: no flap control, 1: steady state flap angle, 2: Proportional flap control}
   # Controller parameters       
-  zeta_pc:            1.0                   # Pitch controller desired damping ratio [-]
-  omega_pc:           0.25                  # Pitch controller desired natural frequency [rad/s]
+  U_pc:               [14,20]
+  zeta_pc:            [2.0,1.0]                  # Pitch controller desired damping ratio [-]
+  omega_pc:           [0.15,0.3]                  # Pitch controller desired natural frequency [rad/s]
   zeta_vs:            0.85                   # Torque controller desired damping ratio [-]
   omega_vs:           0.12                 # Torque controller desired natural frequency [rad/s]
   # Optional - these can be defined, but do not need to be

Tune_Cases/NREL5MW.yaml

@@ -42,8 +42,9 @@ controller_params:
   Fl_Mode:            0                     # Floating specific feedback mode {0: no nacelle velocity feedback, 1: nacelle velocity feedback}
   Flp_Mode:           0                     # Flap control mode {0: no flap control, 1: steady state flap angle, 2: Proportional flap control}
   # Controller parameters 
-  zeta_pc:            0.7                   # Pitch controller desired damping ratio [-]
-  omega_pc:           0.5                   # Pitch controller desired natural frequency [rad/s]
+  U_pc:               [12,18]
+  zeta_pc:            [0.7]                   # Pitch controller desired damping ratio [-]
+  omega_pc:           [0.5]                   # Pitch controller desired natural frequency [rad/s]
   zeta_vs:            0.7                   # Torque controller desired damping ratio [-]
   omega_vs:           0.2                   # Torque controller desired natural frequency [rad/s]
   # Only needed if Flp_Mode > 0