Fixed-Pitch Torque Control in Region 3

Examples/29_marine_hydro_fbp.py

@@ -0,0 +1,139 @@
+'''
+----------- 26_marine_hydro ------------------------
+Run openfast with ROSCO and a MHK turbine
+-----------------------------------------------
+
+
+'''
+
+import os
+from rosco.toolbox.ofTools.case_gen.run_FAST import run_FAST_ROSCO
+from rosco.toolbox.ofTools.case_gen import CaseLibrary as cl
+from rosco.toolbox.ofTools.fast_io import output_processing
+from rosco.toolbox import controller as ROSCO_controller
+from rosco.toolbox import turbine as ROSCO_turbine
+from rosco.toolbox.utilities import write_DISCON
+from rosco.toolbox.inputs.validation import load_rosco_yaml
+import matplotlib.pyplot as plt
+
+'''
+Run MHK turbine in OpenFAST with ROSCO torque controller
+
+
+'''
+
+
+#directories
+this_dir            = os.path.dirname(os.path.abspath(__file__))
+rosco_dir           = os.path.dirname(this_dir)
+example_out_dir     = os.path.join(this_dir,'examples_out')
+os.makedirs(example_out_dir,exist_ok=True)
+
+def main():
+
+    # Input yaml and output directory
+    parameter_filename = os.path.join(this_dir,'Tune_Cases/RM1_MHK_FBP.yaml')
+    run_dir = os.path.join(example_out_dir,'29_MHK/0_baseline')
+    os.makedirs(run_dir,exist_ok=True)
+
+    inps = load_rosco_yaml(parameter_filename)
+    path_params         = inps['path_params']
+    turbine_params      = inps['turbine_params']
+    controller_params   = inps['controller_params']
+
+    # Instantiate turbine, controller, and file processing classes
+    turbine         = ROSCO_turbine.Turbine(turbine_params)
+    controller      = ROSCO_controller.Controller(controller_params)
+
+    # Load turbine data from OpenFAST and rotor performance text file
+    cp_filename = os.path.join(this_dir,path_params['rotor_performance_filename'])
+    turbine.load_from_fast(
+        path_params['FAST_InputFile'],
+        os.path.join(this_dir,path_params['FAST_directory']),
+        rot_source='txt',txt_filename= cp_filename
+        )
+
+    # Tune controller 
+    controller.tune_controller(turbine)
+
+    # Write parameter input file
+    param_file = os.path.join(run_dir,'DISCON.IN')
+    write_DISCON(turbine,
+                 controller,
+                 param_file=param_file, 
+                 txt_filename=cp_filename
+    )
+
+    # # Plot operating 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')
+
+
+    # simulation set up
+    r = run_FAST_ROSCO()
+    r.tuning_yaml   = parameter_filename
+    # r.wind_case_fcn = cl.simp_step  # single step wind input
+    r.wind_case_fcn = cl.power_curve
+    r.wind_case_opts    = {
+        'U': [3.5],
+        'TMax': 100.0,
+        }
+    r.case_inputs = {}
+    # r.fst_vt        = reader.fst_vt
+    # r.controller_params = controller_params
+    r.save_dir      = run_dir
+    r.rosco_dir     = rosco_dir
+
+    r.run_FAST()
+
+    op = output_processing.output_processing()
+    fast_out = op.load_fast_out([os.path.join(run_dir,'RM1_MHK_FBP/power_curve/base/RM1_MHK_FBP_0.out')], tmin=0)
+    fig, axs = plt.subplots(4,1)
+    axs[0].plot(fast_out[0]['Time'],fast_out[0]['Wind1VelX'],label='Flow Speed')
+    axs[1].plot(fast_out[0]['Time'],fast_out[0]['GenSpeed'],label='Gen Speed')
+    axs[2].plot(fast_out[0]['Time'],fast_out[0]['GenTq'],label='Gen Torque')
+    axs[3].plot(fast_out[0]['Time'],fast_out[0]['GenPwr'],label='Gen Power')
+
+
+    # op = output_processing.output_processing()
+    # op2 = output_processing.output_processing()
+
+    # md_out = op.load_fast_out([os.path.join(run_dir,'IEA15MW_cable/power_curve/base/IEA15MW_cable_0.MD.Line1.out')], tmin=0)
+    # local_vars = op2.load_fast_out([os.path.join(run_dir,'IEA15MW_cable/power_curve/base/IEA15MW_cable_0.RO.dbg2')], tmin=0)
+
+    # fig, axs = plt.subplots(4,1)
+    # axs[0].plot(local_vars[0]['Time'],local_vars[0]['CC_DesiredL'],label='CC_DesiredL')
+    # axs[1].plot(local_vars[0]['Time'],local_vars[0]['CC_ActuatedL'],label='CC_ActuatedL')
+    # axs[2].plot(local_vars[0]['Time'],local_vars[0]['CC_ActuatedDL'],label='CC_ActuatedDL')
+    # axs[3].plot(md_out[0]['Time'],md_out[0]['Seg20Lst'],label='Seg20Lst')
+
+    # [a.legend() for a in axs]
+    # [a.grid() for a in axs]
+
+    # if False:
+    #     plt.show()
+    # else:
+    #     plt.savefig(os.path.join(example_out_dir,'22_cable_control.png'))
+
+    # # Check that the last segment of line 1 shrinks by 10 m
+    # # np.testing.assert_almost_equal(md_out[0]['Seg20Lst'][-1] - md_out[0]['Seg20Lst'][0], line_ends[0], 2)
+
+
+
+if __name__=="__main__":
+    main()

Examples/ROSCO_walkthrough.ipynb

@@ -1,6 +1,7 @@
 {
  "cells": [
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -61,6 +62,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -133,6 +135,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -153,6 +156,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -197,6 +201,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -250,6 +255,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -287,6 +293,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -362,6 +369,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -423,6 +431,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -593,6 +602,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -604,6 +614,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -686,6 +697,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {
@@ -724,6 +736,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {
     "slideshow": {

Examples/Test_Cases/MHK_RM1/MHK_RM1_Floating_ElastoDyn.dat

@@ -25,9 +25,9 @@ True          PtfmYDOF    - Platform yaw rotation DOF (flag)
 ---------------------- INITIAL CONDITIONS --------------------------------------
           0   OoPDefl     - Initial out-of-plane blade-tip displacement (meters)
           0   IPDefl      - Initial in-plane blade-tip deflection (meters)
-          10   BlPitch(1)  - Blade 1 initial pitch (degrees)
-          10   BlPitch(2)  - Blade 2 initial pitch (degrees)
-          10   BlPitch(3)  - Blade 3 initial pitch (degrees) [unused for 2 blades]
+          0   BlPitch(1)  - Blade 1 initial pitch (degrees)
+          0   BlPitch(2)  - Blade 2 initial pitch (degrees)
+          0   BlPitch(3)  - Blade 3 initial pitch (degrees) [unused for 2 blades]
           0   TeetDefl    - Initial or fixed teeter angle (degrees) [unused for 3 blades]
           0   Azimuth     - Initial azimuth angle for blade 1 (degrees)
       11.50   RotSpeed    - Initial or fixed rotor speed (rpm)

Examples/Tune_Cases/RM1_MHK_FBP.yaml

@@ -0,0 +1,73 @@
+--- #  ---------------------NREL Generic controller tuning input file -------------------
+    # Written for use with ROSCO_Toolbox tuning procedures
+    # Turbine: RM1 Marine Turbine
+# ------------------------------ OpenFAST PATH DEFINITIONS ------------------------------
+path_params:
+  FAST_InputFile: 'MHK_RM1_Floating.fst'                                                                # Name of *.fst file
+  FAST_directory: '../Test_Cases/MHK_RM1'                    # Main OpenFAST model directory, where the *.fst lives 
+  # Optional
+  rotor_performance_filename:  '../Test_Cases/MHK_RM1/MHK_RM1_Cp_Ct_Cq.txt'                                # Filename for rotor performance text file (if it has been generated by ccblade already)
+
+# -------------------------------- TURBINE PARAMETERS ----------------------------------- 
+turbine_params:
+  rotor_inertia:      92169               # Rotor inertia [kg m^2], {Available in Elastodyn .sum file}
+  rated_rotor_speed:  1.204   # 11.5 rpm               # Rated rotor speed [rad/s]
+  v_min:              0.5                   # Cut-in wind speed [m/s] 
+  v_rated:            2.0                  # Rated wind speed [m/s]
+  v_max:              4.0                  # Cut-out wind speed [m/s], -- Does not need to be exact (JUST ASSUME FOR NOW)
+  max_pitch_rate:     0.1745                # Maximum blade pitch rate [rad/s]
+  max_torque_rate:    1500000.              # Maximum torque rate [Nm/s], {~1/4 VS_RtTq/s}
+  rated_power:        500000              # Rated Power [W]
+  bld_edgewise_freq:  60.2831853             # Blade edgewise first natural frequency [rad/s]
+  bld_flapwise_freq:  0.0                   # Blade flapwise first natural frequency [rad/s]
+  reynolds_ref:       8e6
+  # Optional 
+  # TSR_operational:    # None                # Desired below-rated operational tip speed ratio (Cp-maximizing TSR is used if not defined)
+#------------------------------- CONTROLLER PARAMETERS ----------------------------------
+controller_params:
+  # Controller flags 
+  LoggingLevel:       2                     # {0: write no debug files, 1: write standard output .dbg-file, 2: write standard output .dbg-file and complete avrSWAP-array .dbg2-file
+  F_LPFType:          1                     # {1: first-order low-pass filter, 2: second-order low-pass filter}, [rad/s] (currently filters generator speed and pitch control signals)
+  F_NotchType:        1                     # Notch filter on generator speed and/or tower fore-aft motion (for floating) {0: disable, 1: generator speed, 2: tower-top fore-aft motion, 3: generator speed and tower-top fore-aft motion}
+  IPC_ControlMode:    0                     # Turn Individual Pitch Control (IPC) for fatigue load reductions (pitch contribution) {0: off, 1: 1P reductions, 2: 1P+2P reductions}
+  VS_ControlMode:     4                     # Generator torque control mode in above rated conditions {0: constant torque, 1: constant power, 2: TSR tracking PI control}
+  VS_ConstPower:      0                     # Generator torque control mode in above rated conditions {0: constant torque, 1: constant power, 2: TSR tracking PI control}
+  PC_ControlMode:     0                     # Blade pitch control mode {0: No pitch, fix to fine pitch, 1: active PI blade pitch control}
+  Y_ControlMode:      0                     # Yaw control mode {0: no yaw control, 1: yaw rate control, 2: yaw-by-IPC}
+  SS_Mode:            1                     # Setpoint Smoother mode {0: no setpoint smoothing, 1: introduce setpoint smoothing}
+  WE_Mode:            0                     # Wind speed estimator mode {0: One-second low pass filtered hub height wind speed, 1: Immersion and Invariance Estimator (Ortega et al.)}
+  PS_Mode:            0                     # Pitch saturation mode {0: no pitch saturation, 1: peak shaving, 2: Cp-maximizing pitch saturation, 3: peak shaving and Cp-maximizing pitch saturation}
+  SD_Mode:            0                     # Shutdown mode {0: no shutdown procedure, 1: pitch to max pitch at shutdown}
+  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 
+  U_pc:               [2.3,2.5]
+  interp_type:        sigma
+  zeta_pc:            [0.7,0.7]                   # Pitch controller desired damping ratio [-]
+  omega_pc:           [0.9,0.9]                   # Pitch controller desired natural frequency [rad/s]
+  zeta_vs:            1.0                   # Torque controller desired damping ratio [-]
+  omega_vs:           0.8                   # Torque controller desired natural frequency [rad/s]
+  twr_freq:           3.3404                # Tower natural frequency [rad/s]  # 0.4499 (old value) 3.3404(new value)
+  # twr_freq:           0.061009    # 2P
+  ptfm_freq:          0.6613               # Platform natural frequency [rad/s]  (OC4Hywind Parameters, here)   0.2325 (old value) 0.6613879263 (new value)
+  # Optional         
+  ps_percent:         0.80                  # Percent peak shaving  [%, <= 1 ], {default = 80%}
+  sd_maxpit:          0.4363                # Maximum blade pitch angle to initiate shutdown [rad], {default = bld pitch at v_max}
+  Kp_float:           -0.3897
+  max_torque_factor:  5
+  FBP_ref_mode: 1
+  FBP_P_mode:   0
+  FBP_U:        [2.0, 4.0]
+  FBP_P:        [1.0, 1.0]
+  DISCON:
+    F_NumNotchFilts:  0
+    F_NotchFreqs:  [1.0, 2.4200]  # 2P
+    F_NotchBetaNum: [0.0, 0.0]
+    F_NotchBetaDen: [0.25, 0.25]
+    F_GenSpdNotch_N: 0
+    F_GenSpdNotch_Ind: [1,2]
+    F_TwrTopNotch_N: 0
+    F_TwrTopNotch_Ind: [1,2]
+    F_NotchCornerFreq_GS: [2.42] 
+    F_FlHighPassFreq: 1.5
+    VS_PwrFiltF: 0.05

rosco/controller/rosco_registry/rosco_types.yaml

@@ -1037,6 +1037,12 @@ LocalVariables:
     RotSpeedF:                    
       <<: *real 
       description: Filtered LSS (generator) speed [rad/s].
+    VS_RefSpd:
+      <<: *real
+      description: Generator speed set point of torque controller [rad/s]
+    PC_RefSpd:
+      <<: *real
+      description: Generator speed set point of pitch controller [rad/s]
     GenSpeedF:                    
       <<: *real 
       description: Filtered HSS (generator) speed [rad/s].

rosco/controller/src/ControllerBlocks.f90

@@ -34,6 +34,7 @@ SUBROUTINE ComputeVariablesSetpoints(CntrPar, LocalVar, objInst, DebugVar, ErrVa
         TYPE(DebugVariables),       INTENT(INOUT)       :: DebugVar
         TYPE(ErrorVariables),       INTENT(INOUT)       :: ErrVar
 
+        REAL(DbKi)                 :: VS_RefSpdRaw   ! Temporary variable for applying LPF to the reference speed after it is generated
 
         ! ----- Pitch controller speed and power error -----
 
@@ -62,8 +63,22 @@ SUBROUTINE ComputeVariablesSetpoints(CntrPar, LocalVar, objInst, DebugVar, ErrVa
         IF (CntrPar%VS_ControlMode == 2) THEN
             LocalVar%VS_RefSpd_TSR = (CntrPar%VS_TSRopt * LocalVar%We_Vw_F / CntrPar%WE_BladeRadius) * CntrPar%WE_GearboxRatio
         ELSEIF (CntrPar%VS_ControlMode == 3) THEN
-            LocalVar%VS_GenPwrF = LPFilter(LocalVar%VS_GenPwr, LocalVar%DT,CntrPar%VS_PwrFiltF, LocalVar%FP, LocalVar%iStatus, LocalVar%restart, objInst%instLPF) 
-            LocalVar%VS_RefSpd_TSR = (LocalVar%VS_GenPwrF/CntrPar%VS_Rgn2K)**(1./3.) ! Genspeed reference that doesnt depend on wind speed estimate (https://doi.org/10.2172/1259805)
+            LocalVar%VS_GenPwrF = LPFilter(LocalVar%VS_GenPwr, LocalVar%DT, CntrPar%VS_PwrFiltF, LocalVar%FP, LocalVar%iStatus, LocalVar%restart, objInst%instLPF) 
+            LocalVar%VS_RefSpd_TSR = (LocalVar%VS_GenPwrF/CntrPar%VS_Rgn2K)**(1./3.) ! Genspeed reference that doesn't depend on wind speed estimate (https://doi.org/10.2172/1259805)
+        ELSEIF (CntrPar%VS_ControlMode == 4) THEN
+            IF (CntrPar%FBP_RefMode == 0) THEN
+                ! Use WSE to look up speed reference
+                VS_RefSpdRaw = interp1d(CntrPar%FBP_U, CntrPar%FBP_Omega, LocalVar%WE_Vw, ErrVar)
+            ELSEIF (CntrPar%FBP_RefMode == 1) THEN
+                ! Use LocalVar%GenTq or LocalVar%GenTqMeas, Omega must be expressed as a function of Tau
+                VS_RefSpdRaw = interp1d(CntrPar%FBP_Tau, CntrPar%FBP_Omega, LocalVar%GenTq, ErrVar)
+            ENDIF
+
+            ! VS_RefSpdRaw = (LocalVar%GenTq/CntrPar%VS_Rgn2K)**(1./2.) ! WSE-independent below-rated generator speed reference based on torque instead of power
+            ! VS_RefSpdRaw = min(VS_RefSpdRaw, (CntrPar%VS_RtPwr/(CntrPar%VS_GenEff/100.0))/LocalVar%GenTq) ! Above-rated underspeed reference
+            ! VS_RefSpdRaw = min(VS_RefSpdRaw, CntrPar%VS_RefSpd) ! Saturate to rated speed
+
+            LocalVar%VS_RefSpd_TSR = LPFilter(VS_RefSpdRaw, LocalVar%DT, CntrPar%VS_PwrFiltF, LocalVar%FP, LocalVar%iStatus, LocalVar%restart, objInst%instLPF)
         ELSE
             LocalVar%VS_RefSpd_TSR = CntrPar%VS_RefSpd
         ENDIF 
@@ -97,7 +112,7 @@ SUBROUTINE ComputeVariablesSetpoints(CntrPar, LocalVar, objInst, DebugVar, ErrVa
         LocalVar%VS_RefSpd = max(LocalVar%VS_RefSpd, CntrPar%VS_MinOmSpd)
 
         ! Reference error
-        IF ((CntrPar%VS_ControlMode == 2) .OR. (CntrPar%VS_ControlMode == 3)) THEN
+        IF ((CntrPar%VS_ControlMode == 2) .OR. (CntrPar%VS_ControlMode == 3) .OR. (CntrPar%VS_ControlMode == 4)) THEN
             LocalVar%VS_SpdErr = LocalVar%VS_RefSpd - LocalVar%GenSpeedF
         ENDIF