diff --git a/CMakeLists.txt b/CMakeLists.txt index 5b68871be160788a0063d6bb2fa691c67a286036..dbf159b5568de0cfa372e59f265e2c8d4d579aca 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,6 +11,7 @@ set(DTUWEC_SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src/dtu_we_controller) set(BLADED_SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src/dtu_we_controller_bladed) set(FLAP_SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src/flap_controller_individual_aep_u_f) set(FLAPCYC_SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src/flap_controller_cyclic) +set(YAW_SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src/yaw_controller) set(BUILD_TYPE SHARED) set(LIB static) @@ -28,4 +29,7 @@ message("Configuring flap controller source files: " ${FLAP_SRC_DIR}) add_subdirectory(${FLAP_SRC_DIR}) message("Configuring cyclic flap controller source files: " ${FLAPCYC_SRC_DIR}) -add_subdirectory(${FLAPCYC_SRC_DIR}) \ No newline at end of file +add_subdirectory(${FLAPCYC_SRC_DIR}) + +message("Configuring yaw controller source files: " ${YAW_SRC_DIR}) +add_subdirectory(${YAW_SRC_DIR}) \ No newline at end of file diff --git a/src/yaw_controller/CMakeLists.txt b/src/yaw_controller/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..c0a4db7cbbe3f0f0f5d217554e211421e4666a8d --- /dev/null +++ b/src/yaw_controller/CMakeLists.txt @@ -0,0 +1,25 @@ +# Set the project name +project(yaw_controller LANGUAGES Fortran) + +message("Configuring Sub-project: " ${PROJECT_NAME}) +message("Root Directory is: " ${ROOT_DIR}) + +# set source code +set(SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}) +set(YAWSRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}) + +file(GLOB_RECURSE MODSRC "${ROOT_DIR}/src/modules/global_constants.f90" + "${ROOT_DIR}/src/modules/global_variables.f90" + "${ROOT_DIR}/src/modules/misc_mod.f90" + "${ROOT_DIR}/src/modules/user_defined_types.f90" + ) + +set(SRC ${YAWSRC_DIR}/yaw_controller.f90 + ) +list(APPEND SRC ${MODSRC}) + +message("Including the utils: " "${ROOT_DIR}/utils/cmake/CMakeLists.txt") +include(${ROOT_DIR}/utils/cmake/CMakeLists.txt) + + + diff --git a/src/yaw_controller/file_info.h b/src/yaw_controller/file_info.h new file mode 100644 index 0000000000000000000000000000000000000000..d6046851e6130b92519cab3431d4a694c9959825 --- /dev/null +++ b/src/yaw_controller/file_info.h @@ -0,0 +1,2 @@ +#define FILE_NAME "yaw_controller.dll" +#define FILE_DESCRIPTION "DTU Wind Energy Controller - yaw controller dll" diff --git a/src/yaw_controller/yaw_controller.f90 b/src/yaw_controller/yaw_controller.f90 new file mode 100644 index 0000000000000000000000000000000000000000..6251c81d3b83c353e2d500e1770de7fbcff78bee --- /dev/null +++ b/src/yaw_controller/yaw_controller.f90 @@ -0,0 +1,227 @@ +module YawController +! +! Implements a basic yaw controller. + + + use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan, ieee_is_finite + + +!**************************************************************************************** +! Module variables. +!**************************************************************************************** + +! 32 bit integer. + integer, parameter :: ki_9 = selected_int_kind(9) + +! 64 bit real. + integer, parameter :: kr_dp = selected_real_kind(15, 307) + +! Pi. + real(kind=kr_dp), parameter :: pi = 4.0_kr_dp * atan(1.0_kr_dp) + +! Degrees to radians conversion factor. + real(kind=kr_dp), parameter :: deg2rad = pi / 180.0_kr_dp + +! Radians to degrees conversion factor. + real(kind=kr_dp), parameter :: rad2deg = 180.0_kr_dp / pi + +! Not a Number. +! Defined within init_yaw subroutine. + real(kind=kr_dp) nan + + + type yaw_str + + ! Deltat dll [s]. + real(kind=kr_dp) deltat + + ! Cutoff frequency of yaw error first order low pass filter [1/s]. + real(kind=kr_dp) cutoff_freq + + ! *** VARIABLES *** + ! weight factor of exponential moving average/variance filter + real(kr_dp) :: alpha + + ! Last values of the wind speed and direction. + real(kr_dp) :: wind_speed_ave, previous_wind_speed_ave + real(kr_dp) :: wind_direction_ave, previous_wind_direction_ave + + + ! First iteration flag + logical :: first_iteration = .true. + + end type yaw_str + +! Exchange data from init to update and keep memory. + type(yaw_str) yawst + +contains + + + subroutine init_yaw(array1, array2) bind(C,name='init_yaw') +!DEC$ ATTRIBUTES DLLEXPORT :: init_yaw + + implicit none + + +! Data exchanged with HAWC2. + real(kind=kr_dp), intent(in) :: array1(100) + real(kind=kr_dp), intent(out) :: array2(1) + +! array1 contains: +! 1: Deltat dll [s]. +! 2: Cutoff frequency of yaw error first order low pass filter [1/s]. + + + +! Cutoff frequency converted to angular frequency. + real(kind=kr_dp) :: cutoff_angular_freq + +! Copy parameters from init block to yawst variable. + yawst%deltat = array1( 1) + yawst%cutoff_freq = array1( 2) + + + +! Compute angular frequency of cutoff frequency. + cutoff_angular_freq = 2 * yawst%cutoff_freq * yawst%deltat * pi + +! Computes the exponential filter weighting factor analytically. + yawst%alpha = cos(cutoff_angular_freq) - 1 + sqrt(cos(cutoff_angular_freq)**2 - & + 4*cos(cutoff_angular_freq) + 3) + print*, yawst%alpha +! Define NaN. +! Cannot be a parameter. See https://community.intel.com/t5/Intel-Fortran-Compiler/Set-Parameter-To-NaN/td-p/1132223 + nan = ieee_value(1.0_kr_dp, ieee_quiet_nan) + +! Initialize previous states of filter. + yawst%previous_wind_speed_ave = nan !0.0_kr_dp ! nan + yawst%previous_wind_direction_ave = nan !0.0_kr_dp ! nan + +! Set the state. + +! No output needed for init. + array2 = 0.0_kr_dp + + + + end subroutine init_yaw + + + subroutine update_yaw(array1, array2) bind(C,name='update_yaw') +!DEC$ ATTRIBUTES DLLEXPORT :: update_yaw + + implicit none + + +! Data exchanged with HAWC2. + real(kind=kr_dp), dimension(100), intent(in) :: array1 ! outvec + real(kind=kr_dp), dimension(100), intent(out) :: array2 ! inpvec -> actions + +! array1 contains: +! 1: Time [s] +! 2, 3: Wind speed and direction in non-rotating rotor coordinates [m/s, deg] +! 4: angle of yaw bearing node [rad] +! 5: Yaw setpoint [deg] + +! array2 contains: +! 1: Demanded yaw bearing angle [rad] +! 2: Yaw error [rad] +! 3: Filtered wind speed [m/s] +! 4: Filtered wind direction [rad] +! 5: Instantaneous wind speed [m/s] +! 6: Instantaneous wind direction [rad] +! 7: angle of yaw bearing node [rad] +! 8: yaw setpoint [rad] + +! *** Variables from HAWC2 output block. *** +! general time; Time [s] 1 + real(kind=kr_dp) time + +! wind free_wind_hor_center_pos0 2; +! Instantaneous horizontal wind speed [m/s]. + real(kind=kr_dp) wind_speed_hor_instantaneous + +! Instantaneous wind direction [rad] (converted from deg). + real(kind=kr_dp) wind_direction_instantaneous + +! Current nacelle angle [rad]. +! constraint bearing2 yaw_rot 1 only 1; Rotation of yaw bearing [rad] + real(kind=kr_dp) current_nacelle_angle + +! hawc2farm yaw_setpoint ; Yaw offset setpoint [rad] (converted from deg). + real(kind=kr_dp) yaw_setpoint + + +! *** Return values *** + +! error in yaw alignment [rad] + real(kind=kr_dp) error + +! Demanded nacelle yaw angle [rad] + real(kind=kr_dp) demanded_yaw_bearing_angle + + +! *** Controller code. *** + +! Copy values from the input array. + time = array1(1) + wind_speed_hor_instantaneous = array1(2) + wind_direction_instantaneous = array1(3) * deg2rad + current_nacelle_angle = array1(4) + yaw_setpoint = array1(5) * deg2rad + + + + +! Push most recent filtered values to previous memory slot. + yawst%previous_wind_speed_ave = yawst%wind_speed_ave + yawst%previous_wind_direction_ave = yawst%wind_direction_ave + +! Compute the average wind speed and direction. +! - Since we are using rotor coordinates, yawst%wind_direction_ave is the current yaw error. + if (yawst%first_iteration) then + yawst%wind_speed_ave = wind_speed_hor_instantaneous + yawst%wind_direction_ave = wind_direction_instantaneous + + + yawst%first_iteration = .false. + else + !perform exponential moving average + ! https://stats.stackexchange.com/questions/6874/exponential-weighted-moving-skewness-kurtosis + yawst%wind_speed_ave = (1 - yawst%alpha) * yawst%previous_wind_speed_ave + & + yawst%alpha * wind_speed_hor_instantaneous + yawst%wind_direction_ave = (1 - yawst%alpha) * yawst%previous_wind_direction_ave + & + yawst%alpha * wind_direction_instantaneous + + + endif + + + ! Calculate error between desired and actual yaw alignment. + error = (yawst%wind_direction_ave - yaw_setpoint) + ! jyli: added yaw_bearing_node angle as I believe this compensates any reference frame. + ! jyli: Potential sign error. + ! Calculate yaw command. + demanded_yaw_bearing_angle = current_nacelle_angle - error + + +! Set the action. + array2 = 0.0_kr_dp + array2(1) = demanded_yaw_bearing_angle + +! Set some other outputs of this dll. + array2(2) = error + array2(3) = yawst%wind_speed_ave + array2(4) = yawst%wind_direction_ave + +! Input pass-through. + array2(5) = wind_speed_hor_instantaneous + array2(6) = wind_direction_instantaneous + array2(7) = current_nacelle_angle + array2(8) = yaw_setpoint + + end subroutine update_yaw + + +end module YawController