Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

sigma + ipc #125

Merged
merged 12 commits into from
Apr 25, 2022
13 changes: 12 additions & 1 deletion ROSCO/rosco_registry/rosco_types.yaml
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
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
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
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