Skip to content

Commit

Permalink
sigma + ipc (#125)
Browse files Browse the repository at this point in the history 
* cleanup api change table

* Update inverted notch to move frequency properly

* Saturate inv notch corner frequency at 0

* add sigma function

* Use IPC_Vramp for IPC cut-in

* Add IPC_Vramp DISCON inputs

* update registry

* Update DISCONs

* Update docs for API change

* Fix IPC_Vramp data type

* update comments

Co-authored-by: dzalkind <[email protected]>
  • Loading branch information
@nikhar-abbas @dzalkind
nikhar-abbas and dzalkind authored Apr 25, 2022
1 parent d520192 commit bd0562c
Show file tree
Hide file tree
Showing 14 changed files with 222 additions and 100 deletions.
13 changes: 12 additions & 1 deletion ROSCO/rosco_registry/rosco_types.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,10 @@ ControlParameters:
IPC_ControlMode:
<<: *integer
description: Turn Individual Pitch Control (IPC) for fatigue load reductions (pitch contribution) {0 - off, 1 - 1P reductions, 2 - 1P+2P reductions}
IPC_Vramp:
<<: *real
description: Wind speeds for IPC cut-in sigma function [m/s]
allocatable: True
IPC_IntSat:
<<: *real
description: Integrator saturation (maximum signal amplitude contrbution to pitch from IPC)
Expand Down Expand Up @@ -750,6 +754,14 @@ LocalVariables:
IPC_AxisYaw_2P:
<<: *real
description: Integral of quadrature, 2P
IPC_KI:
<<: *real
description: Integral gain for IPC, after ramp [-]
size: 2
IPC_KP:
<<: *real
description: Proportional gain for IPC, after ramp [-]
size: 2
PC_State:
<<: *integer
description: State of the pitch control system
Expand Down Expand Up @@ -966,7 +978,6 @@ DebugVariables:
<<: *real
description: Yaw component of coleman transformation, 2P


ErrorVariables:
size_avcMSG:
<<: *integer
Expand Down
36 changes: 25 additions & 11 deletions ROSCO/src/Controllers.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ SUBROUTINE PitchControl(avrSWAP, CntrPar, LocalVar, objInst, DebugVar, ErrVar)
DebugVar%PC_PICommand = LocalVar%PC_PitComT
! Find individual pitch control contribution
IF ((CntrPar%IPC_ControlMode >= 1) .OR. (CntrPar%Y_ControlMode == 2)) THEN
CALL IPC(CntrPar, LocalVar, objInst, DebugVar)
CALL IPC(CntrPar, LocalVar, objInst, DebugVar, ErrVar)
ELSE
LocalVar%IPC_PitComF = 0.0 ! THIS IS AN ARRAY!!
END IF
Expand Down Expand Up @@ -105,7 +105,7 @@ SUBROUTINE PitchControl(avrSWAP, CntrPar, LocalVar, objInst, DebugVar, ErrVar)
LocalVar%PitCom(K) = LocalVar%PC_PitComT + LocalVar%FA_PitCom(K)
LocalVar%PitCom(K) = saturate(LocalVar%PitCom(K), LocalVar%PC_MinPit, CntrPar%PC_MaxPit) ! Saturate the command using the pitch satauration limits
LocalVar%PitCom(K) = LocalVar%PC_PitComT + LocalVar%IPC_PitComF(K) ! Add IPC
LocalVar%PitCom(K) = saturate(LocalVar%PitCom(K), LocalVar%PC_MinPit, CntrPar%PC_MaxPit) ! Saturate the command using the absolute pitch angle limits
LocalVar%PitCom(K) = saturate(LocalVar%PitCom(K), CntrPar%PC_MinPit, CntrPar%PC_MaxPit) ! Saturate the command using the absolute pitch angle limits
LocalVar%PitCom(K) = ratelimit(LocalVar%PitCom(K), LocalVar%BlPitch(K), CntrPar%PC_MinRat, CntrPar%PC_MaxRat, LocalVar%DT) ! Saturate the overall command of blade K using the pitch rate limit
END DO

Expand Down Expand Up @@ -281,27 +281,29 @@ SUBROUTINE YawRateControl(avrSWAP, CntrPar, LocalVar, objInst, ErrVar)

END SUBROUTINE YawRateControl
!-------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE IPC(CntrPar, LocalVar, objInst, DebugVar)
SUBROUTINE IPC(CntrPar, LocalVar, objInst, DebugVar, ErrVar)
! Individual pitch control subroutine
! - Calculates the commanded pitch angles for IPC employed for blade fatigue load reductions at 1P and 2P
! - Includes yaw by IPC

USE ROSCO_Types, ONLY : ControlParameters, LocalVariables, ObjectInstances, DebugVariables
USE ROSCO_Types, ONLY : ControlParameters, LocalVariables, ObjectInstances, DebugVariables, ErrorVariables

TYPE(ControlParameters), INTENT(INOUT) :: CntrPar
TYPE(LocalVariables), INTENT(INOUT) :: LocalVar
TYPE(ObjectInstances), INTENT(INOUT) :: objInst
TYPE(DebugVariables), INTENT(INOUT) :: DebugVar
TYPE(DebugVariables), INTENT(INOUT) :: DebugVar
TYPE(ErrorVariables), INTENT(INOUT) :: ErrVar

! Local variables
REAL(DbKi) :: PitComIPC(3), PitComIPCF(3), PitComIPC_1P(3), PitComIPC_2P(3)
INTEGER(IntKi) :: K ! Integer used to loop through turbine blades
INTEGER(IntKi) :: i, K ! Integer used to loop through gains and turbine blades
REAL(DbKi) :: axisTilt_1P, axisYaw_1P, axisYawF_1P ! Direct axis and quadrature axis outputted by Coleman transform, 1P
REAL(DbKi) :: axisTilt_2P, axisYaw_2P, axisYawF_2P ! Direct axis and quadrature axis outputted by Coleman transform, 1P
REAL(DbKi) :: axisYawIPC_1P ! IPC contribution with yaw-by-IPC component
REAL(DbKi) :: Y_MErrF, Y_MErrF_IPC ! Unfiltered and filtered yaw alignment error [rad]


CHARACTER(*), PARAMETER :: RoutineName = 'IPC'

! Body
! Pass rootMOOPs through the Coleman transform to get the tilt and yaw moment axis
CALL ColemanTransform(LocalVar%rootMOOPF, LocalVar%Azimuth, NP_1, axisTilt_1P, axisYaw_1P)
Expand All @@ -316,15 +318,21 @@ SUBROUTINE IPC(CntrPar, LocalVar, objInst, DebugVar)
Y_MErrF = 0.0
Y_MErrF_IPC = 0.0
END IF

! Soft cutin with sigma function
DO i = 1,2
LocalVar%IPC_KP(i) = sigma(LocalVar%WE_Vw, CntrPar%IPC_Vramp(1), CntrPar%IPC_Vramp(2), 0.0_DbKi, CntrPar%IPC_KP(i), ErrVar)
LocalVar%IPC_KI(i) = sigma(LocalVar%WE_Vw, CntrPar%IPC_Vramp(1), CntrPar%IPC_Vramp(2), 0.0_DbKi, CntrPar%IPC_KI(i), ErrVar)
END DO

! Integrate the signal and multiply with the IPC gain
IF ((CntrPar%IPC_ControlMode >= 1) .AND. (CntrPar%Y_ControlMode /= 2)) THEN
LocalVar%IPC_axisTilt_1P = PIController(axisTilt_1P, CntrPar%IPC_KP(1), CntrPar%IPC_KI(1), -CntrPar%IPC_IntSat, CntrPar%IPC_IntSat, LocalVar%DT, 0.0_DbKi, LocalVar%piP, LocalVar%restart, objInst%instPI)
LocalVar%IPC_axisYaw_1P = PIController(axisYawF_1P, CntrPar%IPC_KP(1), CntrPar%IPC_KI(1), -CntrPar%IPC_IntSat, CntrPar%IPC_IntSat, LocalVar%DT, 0.0_DbKi, LocalVar%piP, LocalVar%restart, objInst%instPI)
LocalVar%IPC_axisTilt_1P = PIController(axisTilt_1P, LocalVar%IPC_KP(1), LocalVar%IPC_KI(1), -CntrPar%IPC_IntSat, CntrPar%IPC_IntSat, LocalVar%DT, 0.0_DbKi, LocalVar%piP, LocalVar%restart, objInst%instPI)
LocalVar%IPC_axisYaw_1P = PIController(axisYawF_1P, LocalVar%IPC_KP(1), LocalVar%IPC_KI(1), -CntrPar%IPC_IntSat, CntrPar%IPC_IntSat, LocalVar%DT, 0.0_DbKi, LocalVar%piP, LocalVar%restart, objInst%instPI)

IF (CntrPar%IPC_ControlMode >= 2) THEN
LocalVar%IPC_axisTilt_2P = PIController(axisTilt_2P, CntrPar%IPC_KP(2), CntrPar%IPC_KI(2), -CntrPar%IPC_IntSat, CntrPar%IPC_IntSat, LocalVar%DT, 0.0_DbKi, LocalVar%piP, LocalVar%restart, objInst%instPI)
LocalVar%IPC_axisYaw_2P = PIController(axisYawF_2P, CntrPar%IPC_KP(2), CntrPar%IPC_KI(2), -CntrPar%IPC_IntSat, CntrPar%IPC_IntSat, LocalVar%DT, 0.0_DbKi, LocalVar%piP, LocalVar%restart, objInst%instPI)
LocalVar%IPC_axisTilt_2P = PIController(axisTilt_2P, LocalVar%IPC_KP(2), LocalVar%IPC_KI(2), -CntrPar%IPC_IntSat, CntrPar%IPC_IntSat, LocalVar%DT, 0.0_DbKi, LocalVar%piP, LocalVar%restart, objInst%instPI)
LocalVar%IPC_axisYaw_2P = PIController(axisYawF_2P, LocalVar%IPC_KP(2), LocalVar%IPC_KI(2), -CntrPar%IPC_IntSat, CntrPar%IPC_IntSat, LocalVar%DT, 0.0_DbKi, LocalVar%piP, LocalVar%restart, objInst%instPI)
END IF
ELSE
LocalVar%IPC_axisTilt_1P = 0.0
Expand Down Expand Up @@ -361,6 +369,12 @@ SUBROUTINE IPC(CntrPar, LocalVar, objInst, DebugVar)
DebugVar%axisYaw_2P = axisYaw_2P



! Add RoutineName to error message
IF (ErrVar%aviFAIL < 0) THEN
ErrVar%ErrMsg = RoutineName//':'//TRIM(ErrVar%ErrMsg)
ENDIF

END SUBROUTINE IPC
!-------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE ForeAftDamping(CntrPar, LocalVar, objInst)
Expand Down
46 changes: 32 additions & 14 deletions ROSCO/src/Filters.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,30 +137,47 @@ REAL(DbKi) FUNCTION HPFilter( InputSignal, DT, CornerFreq, FP, iStatus, reset, i

END FUNCTION HPFilter
!-------------------------------------------------------------------------------------------------------------------------------
REAL(DbKi) FUNCTION NotchFilterSlopes(InputSignal, DT, CornerFreq, Damp, FP, iStatus, reset, inst)
REAL(DbKi) FUNCTION NotchFilterSlopes(InputSignal, DT, CornerFreq, Damp, FP, iStatus, reset, inst, Moving)
! Discrete time inverted Notch Filter with descending slopes, G = CornerFreq*s/(Damp*s^2+CornerFreq*s+Damp*CornerFreq^2)
USE ROSCO_Types, ONLY : FilterParameters
TYPE(FilterParameters), INTENT(INOUT) :: FP

REAL(DbKi), INTENT(IN) :: InputSignal
REAL(DbKi), INTENT(IN) :: DT ! time step [s]
REAL(DbKi), INTENT(IN) :: CornerFreq ! corner frequency [rad/s]
REAL(DbKi), INTENT(IN) :: Damp ! Dampening constant
INTEGER(IntKi), INTENT(IN) :: iStatus ! A status flag set by the simulation as follows: 0 if this is the first call, 1 for all subsequent time steps, -1 if this is the final call at the end of the simulation.
INTEGER(IntKi), INTENT(INOUT) :: inst ! Instance number. Every instance of this function needs to have an unique instance number to ensure instances don't influence each other.
LOGICAL(4), INTENT(IN) :: reset ! Reset the filter to the input signal

REAL(DbKi), INTENT(IN) :: InputSignal
REAL(DbKi), INTENT(IN) :: DT ! time step [s]
REAL(DbKi), INTENT(IN) :: CornerFreq ! corner frequency [rad/s]
REAL(DbKi), INTENT(IN) :: Damp ! Dampening constant
INTEGER(IntKi), INTENT(IN) :: iStatus ! A status flag set by the simulation as follows: 0 if this is the first call, 1 for all subsequent time steps, -1 if this is the final call at the end of the simulation.
INTEGER(IntKi), INTENT(INOUT) :: inst ! Instance number. Every instance of this function needs to have an unique instance number to ensure instances don't influence each other.
LOGICAL(4), INTENT(IN) :: reset ! Reset the filter to the input signal
LOGICAL, OPTIONAL, INTENT(IN) :: Moving ! Moving CornerFreq flag

LOGICAL :: Moving_ ! Local version
REAL(DbKi) :: CornerFreq_ ! Local version

Moving_ = .FALSE.
IF (PRESENT(Moving)) Moving_ = Moving

! Saturate Corner Freq at 0
IF (CornerFreq < 0) THEN
CornerFreq_ = 0
ELSE
CornerFreq_ = CornerFreq
ENDIF

! Initialization
IF ((iStatus == 0) .OR. reset) THEN
FP%nfs_OutputSignalLast1(inst) = InputSignal
FP%nfs_OutputSignalLast2(inst) = InputSignal
FP%nfs_InputSignalLast1(inst) = InputSignal
FP%nfs_InputSignalLast2(inst) = InputSignal
FP%nfs_b2(inst) = 2.0 * DT * CornerFreq
ENDIF

IF ((iStatus == 0) .OR. reset .OR. Moving_) THEN
FP%nfs_b2(inst) = 2.0 * DT * CornerFreq_
FP%nfs_b0(inst) = -FP%nfs_b2(inst)
FP%nfs_a2(inst) = Damp*DT**2.0*CornerFreq**2.0 + 2.0*DT*CornerFreq + 4.0*Damp
FP%nfs_a1(inst) = 2.0*Damp*DT**2.0*CornerFreq**2.0 - 8.0*Damp
FP%nfs_a0(inst) = Damp*DT**2.0*CornerFreq**2.0 - 2*DT*CornerFreq + 4.0*Damp
FP%nfs_a2(inst) = Damp*DT**2.0*CornerFreq_**2.0 + 2.0*DT*CornerFreq_ + 4.0*Damp
FP%nfs_a1(inst) = 2.0*Damp*DT**2.0*CornerFreq_**2.0 - 8.0*Damp
FP%nfs_a0(inst) = Damp*DT**2.0*CornerFreq_**2.0 - 2*DT*CornerFreq_ + 4.0*Damp
ENDIF

NotchFilterSlopes = 1.0/FP%nfs_a2(inst) * (FP%nfs_b2(inst)*InputSignal + FP%nfs_b0(inst)*FP%nfs_InputSignalLast1(inst) &
Expand Down Expand Up @@ -270,7 +287,8 @@ SUBROUTINE PreFilterMeasuredSignals(CntrPar, LocalVar, DebugVar, objInst, ErrVar
! Blade root bending moment for IPC
DO K = 1,LocalVar%NumBl
IF ((CntrPar%IPC_ControlMode > 0) .OR. (CntrPar%Flp_Mode == 3)) THEN
LocalVar%RootMOOPF(K) = NotchFilterSlopes(LocalVar%rootMOOP(K), LocalVar%DT, LocalVar%RotSpeedF, 0.7_DbKi, LocalVar%FP, LocalVar%iStatus, LocalVar%restart, objInst%instNotchSlopes)
! Moving inverted notch at rotor speed to isolate 1P
LocalVar%RootMOOPF(K) = NotchFilterSlopes(LocalVar%rootMOOP(K), LocalVar%DT, LocalVar%RotSpeedF, 0.7_DbKi, LocalVar%FP, LocalVar%iStatus, LocalVar%restart, objInst%instNotchSlopes, .TRUE.)
ELSEIF ( CntrPar%Flp_Mode == 2 ) THEN
! Filter Blade root bending moments
LocalVar%RootMOOPF(K) = SecLPFilter(LocalVar%rootMOOP(K),LocalVar%DT, CntrPar%F_FlpCornerFreq(1), CntrPar%F_FlpCornerFreq(2), LocalVar%FP, LocalVar%iStatus, LocalVar%restart,objInst%instSecLPF)
Expand Down
37 changes: 37 additions & 0 deletions ROSCO/src/Functions.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -510,6 +510,43 @@ REAL(DbKi) FUNCTION AeroDynTorque(RotSpeed, BldPitch, LocalVar, CntrPar, PerfDat

END FUNCTION AeroDynTorque

!-------------------------------------------------------------------------------------------------------------------------------
REAL(DbKi) FUNCTION sigma(x, x0, x1, y0, y1, ErrVar)
! Generic sigma function
USE ROSCO_Types, ONLY : ErrorVariables
IMPLICIT NONE

! Inputs
TYPE(ErrorVariables), INTENT(INOUT) :: ErrVar

REAL(DbKi), Intent(IN) :: x, x0, x1
REAL(DbKi), Intent(IN) :: y0, y1

! Local
REAL(DbKi) :: a3, a2, a1, a0

CHARACTER(*), PARAMETER :: RoutineName = 'sigma'

a3 = 2/(x0-x1)**3
a2 = -3*(x0+x1)/(x0-x1)**3
a1 = 6*x1*x0/(x0-x1)**3
a0 = (x0-3*x1)*x0**2/(x0-x1)**3

IF (x < x0) THEN
sigma = y0
ELSEIF (x > x1) THEN
sigma = y1
ELSE
sigma = (a3*x**3 + a2*x**2 + a1*x + a0)*(y1-y0) + y0
ENDIF

! Add RoutineName to error message
IF (ErrVar%aviFAIL < 0) THEN
ErrVar%ErrMsg = RoutineName//':'//TRIM(ErrVar%ErrMsg)
ENDIF

END FUNCTION sigma


!-------------------------------------------------------------------------------------------------------------------------------
! Copied from NWTC_IO.f90
Expand Down
Loading

0 comments on commit bd0562c

Please sign in to comment.