diff --git a/DISCON/DISCON_gwin32.dll b/DISCON/DISCON_gwin32.dll index 5e07fb1e..837fe669 100644 Binary files a/DISCON/DISCON_gwin32.dll and b/DISCON/DISCON_gwin32.dll differ diff --git a/Parameter_files/DTU10MW/DISCON.IN b/Parameter_files/DTU10MW/DISCON.IN index b2d7f3ae..ea122128 100644 --- a/Parameter_files/DTU10MW/DISCON.IN +++ b/Parameter_files/DTU10MW/DISCON.IN @@ -1,10 +1,10 @@ -0 = User Variable 1: Yaw control mode (0 = no yaw control, 1 = yaw rate control, 2 = yaw-by-IPC) +0 = User Variable 1: Yaw control mode (0 = no yaw control, 1 = yaw rate control, 2 = yaw-by-IPC) 0.005235988 = User Variable 2: Yaw rate, [rad/s] 26.179938 = User Variable 3: VS_CtInSp, cut-in speed, [rad/s] 150000.0 = User Variable 4: VS_MaxRat, maximum torque change rate, [Nm/s] 50.26548 = User Variable 5: PC_RefSpd, generator speed setpoint for the pitch controller, [rad/s] -0 = User Variable 6: IPC mode for load reductions (0 = off, 1 = on) -0 = User Variable 7: Above rated control (0 = constant torque control, 1 = constant power control) -0.0 = User Variable 8: +0 = User Variable 6: IPC mode for load reductions (0 = off, 1 = on) +0 = User Variable 7: Above rated control (0 = constant torque control, 1 = constant power control) +0.00000 = User Variable 8: Yaw alignment error, setpoint [rad] 0.0 = User Variable 9: -0.0 = User Variable 10: \ No newline at end of file +0.0 = User Variable 10: diff --git a/Parameter_files/DTU10MW/PitchGains.IN b/Parameter_files/DTU10MW/PitchGains.IN index be49fa4b..e381252f 100644 --- a/Parameter_files/DTU10MW/PitchGains.IN +++ b/Parameter_files/DTU10MW/PitchGains.IN @@ -1,9 +1,12 @@ -12 -0.0719 0.1255 0.1638 0.1958 0.2752 0.2985 0.3206 0.3417 0.3622 0.3819 0.4011 0.4196 --0.0338 -0.0172 -0.0141 -0.0118 -0.0079 -0.0066 -0.0055 -0.005 -0.0035 -0.0031 -0.0026 -0.0017 --0.0044 -0.0039 -0.0018 -0.0016 -0.0016 -0.0017 -0.0018 -0.0019 -0.003 -0.003 -0.003 -0.0032 -0 0 0 0 0 0 0 0 0 0 0 0 -0 0 0 0 0 0 0 0 0 0 0 0 -1 --3.48e+04 --2.26e+04 \ No newline at end of file +13 ! Pitch control: amount of gain-scheduled parameters +0 0.0349 0.0698 0.1047 0.1396 0.1745 0.2094 0.2443 0.2793 0.3142 0.3491 0.384 0.4189 ! Pitch control: pitch angle breakpoints (rad) +-0.0105 -0.0104 -0.0102 -0.0099 -0.0095 -0.0091 -0.0086 -0.0081 -0.0076 -0.0071 -0.0066 -0.0061 -0.0056 ! Pitch control: proportional gain (s) +-0.0028 -0.0028 -0.0027 -0.0027 -0.0026 -0.0024 -0.0023 -0.0022 -0.002 -0.0019 -0.0018 -0.0016 -0.0015 ! Pitch control: integral gain (-) +0 0 0 0 0 0 0 0 0 0 0 0 0 ! Pitch control: derivative gain (1/s) +0 0 0 0 0 0 0 0 0 0 0 0 0 ! Pitch control: first-order filter coefficient for derivative action +1 ! Generator torque control: amount of gain-scheduled parameters +-27338.24 ! Generator torque control: proportional gain (s/rad Nm) +-6134.68 ! Generator torque control: integral gain (1/rad Nm) +1 ! Yaw-by-IPC control: amount of gain-scheduled parameters +-0.16 ! Yaw-by-IPC control: proportional gain (1/rad Nm) +-0.002 ! Yaw-by-IPC control: integral gain (1/(rad s) Nm) diff --git a/Parameter_files/DTU10MW/ServoDynBladedInterface.dat b/Parameter_files/DTU10MW/ServoDynBladedInterface.dat index d9d2622a..e8d84ef0 100644 --- a/Parameter_files/DTU10MW/ServoDynBladedInterface.dat +++ b/Parameter_files/DTU10MW/ServoDynBladedInterface.dat @@ -1,10 +1,10 @@ ----------------------- BLADED INTERFACE ---------------------------------------- [used only with Bladed Interface] -"../Controller/DISCON_gwin32.dll" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] -"DISCON.IN" DLL_InFile - Name of input file sent to the DLL (-) [used only with Bladed Interface] +---------------------- BLADED INTERFACE ---------------------------------------- +".\ServoData\DISCON_gwin32.dll" DLL_FileName - Name/location of the dynamic library {.dll [Windows] or .so [Linux]} in the Bladed-DLL format (-) [used only with Bladed Interface] +"DISCON.IN" DLL_InFile - Name of input file sent to the DLL [used only with Bladed Interface] (-) "DISCON" DLL_ProcName - Name of procedure in DLL to be called (-) [case sensitive; used only with DLL Interface] -"default" DLL_DT - Communication interval for dynamic library (s) (or "default") [used only with Bladed Interface] -false DLL_Ramp - Whether a linear ramp should be used between DLL_DT time steps [introduces time shift when true] (flag) [used only with Bladed Interface] - 9999.9 BPCutoff - Cuttoff frequency for low-pass filter on blade pitch from DLL (Hz) [used only with Bladed Interface] +"default" DLL_DT - Communication interval for dynamic library [used only with Bladed Interface] (s) (or "default") +false DLL_Ramp - Whether a linear ramp should be used between DLL_DT time steps [introduces time shift when true] (flag) +9999.9 BPCutoff - Cuttoff frequency for low-pass filter on blade pitch from DLL (Hz) 0 NacYaw_North - Reference yaw angle of the nacelle when the upwind end points due North (deg) [used only with Bladed Interface] 1 Ptch_Cntrl - Record 28: Use individual pitch control {0: collective pitch; 1: individual pitch control} (switch) [used only with Bladed Interface] 0 Ptch_SetPnt - Record 5: Below-rated pitch angle set-point (deg) [used only with Bladed Interface] @@ -12,7 +12,7 @@ false DLL_Ramp - Whether a linear ramp should be used between DLL_DT 90.0 Ptch_Max - Record 7: Maximum pitch angle (deg) [used only with Bladed Interface] -10.0 PtchRate_Min - Record 8: Minimum pitch rate (most negative value allowed) (deg/s) [used only with Bladed Interface] 10.0 PtchRate_Max - Record 9: Maximum pitch rate (deg/s) [used only with Bladed Interface] - 80.1 Gain_OM - Record 16: Optimal mode gain (Nm/(rad/s)^2) [used only with Bladed Interface] + 104.1 Gain_OM - Record 16: Optimal mode gain (Nm/(rad/s)^2) [used only with Bladed Interface] 300.0 GenSpd_MinOM - Record 17: Minimum generator speed (rpm) [used only with Bladed Interface] 405.0 GenSpd_MaxOM - Record 18: Optimal mode maximum speed (rpm) [used only with Bladed Interface] 405.0 GenSpd_Dem - Record 19: Demanded generator speed above rated (rpm) [used only with Bladed Interface] diff --git a/Source/DISCON.f90 b/Source/DISCON.f90 index e13e20fa..a648469a 100644 --- a/Source/DISCON.f90 +++ b/Source/DISCON.f90 @@ -40,12 +40,15 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM REAL(4) :: Azimuth ! Rotor azimuth angle [rad]. REAL(4) :: BlPitch(3) ! Current values of the blade pitch angles [rad]. -REAL(4), PARAMETER :: CornerFreq = 1.570796 ! Corner frequency (-3dB point) in the first-order low-pass filter, [rad/s] +REAL(4), PARAMETER :: CornerFreq = 2.5132741 ! Corner frequency (-3dB point) in the first-order low-pass filter, [rad/s] REAL(4) :: DT ! Time step [s]. REAL(4) :: ElapTime ! Elapsed time since the last call to the controller [s]. REAL(4) :: GenSpeed ! Current HSS (generator) speed [rad/s]. REAL(4) :: GenSpeedF ! Filtered HSS (generator) speed [rad/s]. REAL(4) :: GenTrq ! Electrical generator torque, [Nm]. +REAL(4) :: GenTrqAr ! Electrical generator torque, for above-rated PI-control [Nm]. +REAL(4), PARAMETER :: GenTrqArSatMax = 212900 ! Above rated generator torque PI control saturation +REAL(4) :: GenTrqBr ! Electrical generator torque, for below-rated PI-control [Nm]. REAL(4) :: HorWindV ! Horizontal wind speed at hub-height, [m/s]. REAL(4), PARAMETER :: IPC_KI = 8E-10 ! Integral gain for the individual pitch controller, [-]. INTEGER(4), SAVE :: IPC_ControlMode ! Turn Individual Pitch Control (IPC) for fatigue load reductions (pitch contribution) on = 1/off = 1 @@ -61,6 +64,7 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM REAL(4), SAVE :: LastTime ! Last time this DLL was called, [s]. REAL(4), SAVE :: LastTimePC ! Last time the pitch controller was called, [s]. REAL(4), SAVE :: LastTimeVS ! Last time the torque controller was called, [s]. +INTEGER(4) :: PC_ControlModeNow = 0 ! Current control mode REAL(4) :: PC_GK ! Current value of the gain correction factor, used in the gain scheduling law of the pitch controller, [-]. INTEGER(4), SAVE :: PC_GS_n ! Amount of gain-scheduling table entries REAL(4), DIMENSION(:), ALLOCATABLE, SAVE :: PC_GS_angles ! Gain-schedule table: pitch angles @@ -77,48 +81,58 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM REAL(4) :: PC_MaxRat ! Maximum pitch rate (in absolute value) in pitch controller, [rad/s]. REAL(4) :: PC_MinPit ! Minimum physical pitch limit, [rad]. REAL(4) :: PC_MinRat ! Minimum pitch rate (in absolute value) in pitch controller, [rad/s]. +REAL(4) :: PC_PwrErr ! Power error with respect to rated power [W] REAL(4), SAVE :: PC_RefSpd ! Desired (reference) HSS speed for pitch controller, [rad/s]. REAL(4) :: PC_RtTq99 ! 99% of the rated torque value, using for switching between pitch and torque control, [Nm]. REAL(4) :: PC_SetPnt ! Fine pitch angle, [rad]. REAL(4) :: PC_SpdErr ! Current speed error (pitch control) [rad/s]. REAL(4), SAVE :: PitCom(3) ! Commanded pitch of each blade the last time the controller was called, [rad]. REAL(4), SAVE :: PitComT ! Total command pitch based on the sum of the proportional and integral terms, [rad]. +REAL(4) :: PitComTF ! Total collective pitch command, filtered for gain-scheduling REAL(4), SAVE :: PitComT_IPC(3) ! Total command pitch based on the sum of the proportional and integral terms, including IPC term [rad]. REAL(4), PARAMETER :: R2D = 57.295780 ! Factor to convert radians to degrees. REAL(4) :: rootMOOP(3) ! Blade root out of plane bending moments, [Nm]. REAL(4), PARAMETER :: RPS2RPM = 9.5492966 ! Factor to convert radians per second to revolutions per minute. REAL(4) :: Time ! Current simulation time, [s]. INTEGER(4), SAVE :: VS_ControlMode ! Generator torque control mode in above rated conditions, 0 = constant torque / 1 = constant power +INTEGER(4) :: VS_ControlModeNow = 0 ! Current control mode REAL(4), SAVE :: VS_CtInSp ! Transitional generator speed (HSS side) between regions 1 and 1 1/2, [rad/s]. +REAL(4) :: VS_GenPwr ! Measured generator power [W] INTEGER(4), SAVE :: VS_n ! Number of controller gains REAL(4), DIMENSION(:), ALLOCATABLE, SAVE :: VS_KP ! Proportional gain for generator PI torque controller, used in the transitional 2.5 region REAL(4), DIMENSION(:), ALLOCATABLE, SAVE :: VS_KI ! Integral gain for generator PI torque controller, used in the transitional 2.5 region REAL(4) :: VS_MaxOM ! Optimal mode maximum speed, [rad/s]. REAL(4), SAVE :: VS_MaxRat ! Maximum torque rate (in absolute value) in torque controller, [Nm/s]. REAL(4) :: VS_MaxTq ! Maximum generator torque in Region 3 (HSS side), [Nm]. -- chosen to be 10% above VS_RtTq -REAL(4) :: VS_MinOM ! Optimal mode minimum speed, [rad/s]. -REAL(4) :: VS_Rgn2K ! Generator torque constant in Region 2 (HSS side), N-m/(rad/s)^2. +REAL(4) :: VS_MinTq ! Minimum generator (HSS side), [Nm]. +REAL(4) :: VS_MinOM ! Optimal mode minimum speed, [rad/s] +REAL(4) :: VS_Rgn2K ! Generator torque constant in Region 2 (HSS side), N-m/(rad/s)^2 REAL(4), SAVE :: VS_Rgn2MaxTq ! Maximum torque at the end of the below-rated region 2, [Nm] +REAL(4), SAVE :: VS_Rgn2MinTq ! Minimum torque at the beginning of the below-rated region 2, [Nm] REAL(4), SAVE :: VS_Rgn3MP ! Minimum pitch angle at which the torque is computed as if we are in region 3 regardless of the generator speed, [rad]. -- chosen to be 1.0 degree above PC_SetPnt REAL(4), SAVE :: VS_RtPwr ! Wind turbine rated power [W] REAL(4) :: VS_RtTq ! Rated torque, [Nm]. REAL(4) :: VS_RtSpd ! Rated generator speed [rad/s] REAL(4), SAVE :: VS_Slope15 ! Torque/speed slope of region 1 1/2 cut-in torque ramp , [Nm/(rad/s)]. -REAL(4) :: VS_SpdErr ! Current speed error (generator torque control) [rad/s]. +REAL(4) :: VS_SpdErrAr ! Current speed error (generator torque control) [rad/s]. +REAL(4) :: VS_SpdErrBr ! Current speed error (generator torque control) [rad/s]. REAL(4), SAVE :: Y_AccErr ! Accumulated yaw error [rad]. INTEGER(4), SAVE :: Y_ControlMode ! Yaw control mode: (0 = no yaw control, 1 = yaw rate control, 2 = yaw-by-IPC) REAL(4) :: Y_ErrLPFFast ! Filtered yaw error by fast low pass filter [rad]. REAL(4) :: Y_ErrLPFSlow ! Filtered yaw error by slow low pass filter [rad]. -REAL(4), PARAMETER :: Y_ErrThresh = 1.745329252 ! Error threshold [rad]. Turbine begins to yaw when it passes this. (104.71975512). -REAL(4), SAVE :: Y_YawRate ! Yaw rate [rad/s]. -REAL(4) :: Y_MErr ! Measured yaw error [rad]. -REAL(4), PARAMETER :: Y_omegaLPFast = 1.0 ! Corner frequency fast low pass filter, [Hz]. -REAL(4), PARAMETER :: Y_omegaLPSlow = 0.05 ! Corner frequency slow low pass filter, 1/60 [Hz]. -REAL(4), SAVE :: Y_YawEndT ! Yaw end time, [s]. Indicates the time up until which yaw is active with a fixed rate. - -REAL(4), SAVE :: testValue ! TestValue -REAL(4) :: PC_KPtest -REAL(4) :: PC_KItest +REAL(4), PARAMETER :: Y_ErrThresh = 1.745329252 ! Error threshold [rad]. Turbine begins to yaw when it passes this. (104.71975512) +INTEGER(4), SAVE :: Y_IPC_n ! Number of controller gains (yaw-by-IPC) +REAL(4), DIMENSION(:), ALLOCATABLE, SAVE :: Y_IPC_KP ! Yaw-by-IPC proportional controller gain Kp +REAL(4), DIMENSION(:), ALLOCATABLE, SAVE :: Y_IPC_KI ! Yaw-by-IPC integral controller gain Ki +REAL(4), SAVE :: Y_YawRate ! Yaw rate [rad/s] +REAL(4) :: Y_M ! Yaw misalignment with wind direction, measured [rad] +REAL(4) :: Y_MErr ! Measured yaw error, measured [rad] +REAL(4), SAVE :: Y_MErrSet ! Yaw alignment error, setpoint [rad] +REAL(4), PARAMETER :: Y_omegaLPFast = 1.0 ! Corner frequency fast low pass filter, [Hz] +REAL(4), PARAMETER :: Y_omegaLPSlow = 0.05 ! Corner frequency slow low pass filter, 1/60 [Hz] +REAL(4), SAVE :: Y_YawEndT ! Yaw end time, [s]. Indicates the time up until which yaw is active with a fixed rate + +! REAL(4), SAVE :: testValue ! TestValue INTEGER(4) :: I ! Generic index. INTEGER(4) :: 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. @@ -151,13 +165,14 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM PC_MinRat = avrSWAP(8) PC_MaxRat = avrSWAP(9) VS_RtPwr = avrSWAP(13) +VS_GenPwr = avrSWAP(15) VS_Rgn2K = avrSWAP(16) VS_MinOM = avrSWAP(17) VS_MaxOM = avrSWAP(18) VS_RtSpd = avrSWAP(19) GenSpeed = avrSWAP(20) VS_RtTq = avrSWAP(22) -Y_MErr = avrSWAP(24) +Y_M = avrSWAP(24) HorWindV = avrSWAP(27) rootMOOP(1) = avrSWAP(30) rootMOOP(2) = avrSWAP(31) @@ -168,9 +183,11 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM NumBl = NINT(avrSWAP(61)) PC_RtTq99 = VS_RtTq*0.99 +VS_MinTq = 0.0 VS_MaxTq = VS_RtTq*1.1 +VS_Rgn2MinTq = VS_Rgn2K*VS_MinOM**2 VS_Rgn2MaxTq = VS_Rgn2K*VS_MaxOM**2 -VS_Rgn3MP = PC_SetPnt + 1.0/R2D +VS_Rgn3MP = PC_SetPnt + 2.0/R2D ! Convert C character arrays to Fortran strings: @@ -194,9 +211,8 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM ! Inform users that we are using this user-defined routine: aviFAIL = 1 - ErrMsg = 'Running the Delft Research Controller (DRC) for the NREL' // & - '5MW baseline wind turbine from DISCON.dll.' // & - 'Delft Univeristy of Technology' // & + ErrMsg = 'Running the Delft Research Controller (DRC)\n' // & + 'Delft Univeristy of Technology\n' // & 'Based on baseline NREL5MW controller by J. Jonkman' ! Read user defined parameter file @@ -214,10 +230,11 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM PC_RefSpd = avrSWAP(124) IPC_ControlMode = avrSWAP(125) VS_ControlMode = avrSWAP(126) + Y_MErrSet = avrSWAP(127) ! Initialize testValue (debugging variable) - testValue = 0.0 + ! testValue = 0.4 ! Determine some torque control parameters not specified directly: @@ -263,6 +280,14 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM ALLOCATE(VS_KI(VS_n)) READ(UnPitchGains,*) VS_KI + + READ(UnPitchGains, *) Y_IPC_n + + ALLOCATE(Y_IPC_KP(Y_IPC_n)) + READ(UnPitchGains,*) Y_IPC_KP + + ALLOCATE(Y_IPC_KI(Y_IPC_n)) + READ(UnPitchGains,*) Y_IPC_KI !.............................................................................................................................. ! Check validity of input parameters: @@ -397,8 +422,8 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM IF (DbgOut) THEN OPEN (UnDb, FILE=TRIM(RootName)//'.dbg', STATUS='REPLACE') - WRITE (UnDb,'(A)') ' Time ' //Tab//'PitComT ' //Tab//'PC_SpdError ' //Tab//'PC_KP ' //Tab//'PC_KI ' //Tab//'Y_MErr ' //Tab//'rootMOOP(1) '//Tab//'VS_RtPwr '//Tab//'GenTrq ' - WRITE (UnDb,'(A)') ' (sec) ' //Tab//'(rad) ' //Tab//'(rad/s) '//Tab//'(-) ' //Tab//'(-) ' //Tab//'(rad) ' //Tab//'(?) ' //Tab//'(W) '//Tab//'(Nm) ' + WRITE (UnDb,'(A)') ' Time ' //Tab//'PitComT ' //Tab//'PC_SpdError ' //Tab//'PC_KP ' //Tab//'PC_KI ' //Tab//'Y_M ' //Tab//'rootMOOP(1) '//Tab//'VS_RtPwr '//Tab//'GenTrq '//Tab//'VS_ControlModeNow '//Tab//'PC_ControlModeNow ' + WRITE (UnDb,'(A)') ' (sec) ' //Tab//'(rad) ' //Tab//'(rad/s) '//Tab//'(-) ' //Tab//'(-) ' //Tab//'(rad) ' //Tab//'(?) ' //Tab//'(W) '//Tab//'(Nm) '//Tab//'(-) ' //Tab//'(-) ' OPEN(UnDb2, FILE=TRIM(RootName)//'.dbg2', STATUS='REPLACE') WRITE(UnDb2,'(/////)') @@ -437,8 +462,11 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM ! Filter the HSS (generator) speed measurement: ! Apply Low-Pass Filter - GenSpeedF = LPFilter(GenSpeed, DT, CornerFreq, iStatus, .FALSE., 1) ! This is the first instance of LPFilter + GenSpeedF = SecLPFilter(GenSpeed, DT, CornerFreq, 0.7, iStatus, 1) ! This is the first instance of a second order LPFilter + ! Calculate yaw-alignment error + Y_MErr = Y_M + Y_MErrSet + !.............................................................................................................................. ! VARIABLE-SPEED TORQUE CONTROL: !.............................................................................................................................. @@ -449,25 +477,35 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM ! Compute the generator torque, which depends on which region we are in: - VS_SpdErr = VS_RtSpd - GenSpeedF ! Current speed error + VS_SpdErrAr = VS_RtSpd - GenSpeedF ! Current speed error - Above-rated PI-control + VS_SpdErrBr = VS_MinOM - GenSpeedF ! Current speed error - Below-rated PI-control IF (PitComT >= VS_Rgn3MP) THEN ! We are in region 3 + GenTrqAr = PIController(VS_SpdErrAr, VS_KP(1), VS_KI(1), VS_Rgn2MaxTq, GenTrqArSatMax, DT, GenTrqArSatMax, .TRUE., 1) + GenTrqBr = PIController(VS_SpdErrBr, VS_KP(1), VS_KI(1), VS_MinTq, VS_Rgn2MinTq, DT, VS_Rgn2MinTq, .TRUE., 4) IF (VS_ControlMode == 1) THEN ! Constant power tracking GenTrq = VS_RtPwr/GenSpeedF + VS_ControlModeNow = 6 ELSE ! Constant torque tracking GenTrq = VS_RtTq + VS_ControlModeNow = 5 END IF ELSE - GenTrq = PIController(VS_SpdErr, VS_KP(1), VS_KI(1), VS_Rgn2MaxTq, VS_RtTq, DT, VS_Rgn2MaxTq, 1) - IF (GenTrq >= VS_Rgn2MaxTq*1.01) THEN + GenTrqAr = PIController(VS_SpdErrAr, VS_KP(1), VS_KI(1), VS_Rgn2MaxTq, GenTrqArSatMax, DT, VS_Rgn2MaxTq, .FALSE., 1) + GenTrqBr = PIController(VS_SpdErrBr, VS_KP(1), VS_KI(1), VS_MinTq, VS_Rgn2MinTq, DT, VS_Rgn2MinTq, .FALSE., 4) + IF (GenTrqAr >= VS_Rgn2MaxTq*1.01) THEN + GenTrq = GenTrqAr + VS_ControlModeNow = 3 + CONTINUE + ELSEIF (GenTrqBr <= VS_Rgn2MinTq*0.99) THEN ! We are in region 1 1/2 + GenTrq = GenTrqBr + VS_ControlModeNow = 1 CONTINUE - ELSEIF (GenSpeedF <= VS_CtInSp) THEN ! We are in region 1 - torque is zero - GenTrq = 0.0 - ELSEIF (GenSpeedF < VS_MinOM) THEN ! We are in region 1 1/2 - linear ramp in torque from zero to optimal - GenTrq = VS_Slope15*(GenSpeedF - VS_CtInSp) ELSEIF (GenSpeedF < VS_MaxOM) THEN ! We are in region 2 - optimal torque is proportional to the square of the generator speed GenTrq = VS_Rgn2K*GenSpeedF*GenSpeedF + VS_ControlModeNow = 2 ELSE ! We are in region 2 1/2 - simple induction generator transition region GenTrq = VS_Rgn2MaxTq + VS_ControlModeNow = 4 END IF END IF @@ -497,12 +535,15 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM ! Pitch control !.............................................................................................................................. - IF (VS_ControlMode == 0 .AND. GenTrq > PC_RtTq99) THEN + IF (VS_ControlMode == 0 .AND. GenTrq >= PC_RtTq99) THEN PC_MaxPitVar = PC_MaxPit - ELSEIF (VS_ControlMode == 1 .AND. GenTrq >= VS_RtPwr/GenSpeedF*0.99) THEN + PC_ControlModeNow = 2 + ELSEIF (VS_ControlMode == 1 .AND. GenTrqAr >= GenTrqArSatMax*0.99) THEN PC_MaxPitVar = PC_MaxPit + PC_ControlModeNow = 3 ELSE PC_MaxPitVar = PC_SetPnt + PC_ControlModeNow = 1 END IF ! Compute the elapsed time since the last call to the controller: @@ -511,28 +552,31 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM ! Compute the gain scheduling correction factor based on the previously ! commanded pitch angle for blade 1: - + + ! PitComTF = LPFilter(PitComT, DT, 1.0, iStatus, .FALSE., 1) PC_KP = interp1d(PC_GS_angles, PC_GS_kp, PitComT) PC_KI = interp1d(PC_GS_angles, PC_GS_ki, PitComT) PC_KD = interp1d(PC_GS_angles, PC_GS_kd, PitComT) PC_TF = interp1d(PC_GS_angles, PC_GS_tf, PitComT) - PC_KPtest = interp1d(PC_GS_angles, PC_GS_kp, testValue) - PC_KItest = interp1d(PC_GS_angles, PC_GS_ki, testValue) - testValue = testValue + 0.001 ! Compute the current speed error and its integral w.r.t. time; saturate the ! integral term using the pitch angle limits: - PC_SpdErr = PC_RefSpd - GenSpeedF ! Current speed error + PC_SpdErr = PC_RefSpd - GenSpeedF ! Speed error + PC_PwrErr = VS_RtPwr - VS_GenPwr ! Power error ! Compute the pitch commands associated with the proportional and integral ! gains: - - PitComT = PIController(PC_SpdErr, PC_KP, PC_KI, PC_SetPnt, PC_MaxPitVar, ElapTime, PC_SetPnt, 2) + DFController(PC_SpdErr, PC_KD, PC_TF, DT, 1) + + ! PitComT = NotchFilter(PC_SpdErr, DT, 1.59, 0.01, 0.2, iStatus, 1) + PitComT = PIController(PC_SpdErr, PC_KP, PC_KI, PC_SetPnt, PC_MaxPitVar, DT, PC_SetPnt, .FALSE., 2) ! + DFController(PC_SpdErr, PC_KD, PC_TF, DT, 1) + IF (VS_ControlMode == 1) THEN + PitComT = PitComT + PIController(PC_PwrErr, -4.0E-09, -4.0E-09, PC_SetPnt, PC_MaxPitVar, DT, PC_SetPnt, .FALSE., 5) + END IF ! Individual pitch control - IF (IPC_ControlMode == 1) THEN - CALL IPC(rootMOOP, Azimuth, IPC_phi, Y_MErr, DT, IPC_KI, IPC_omegaHP, IPC_omegaLP, IPC_omegaNotch, IPC_zetaHP, IPC_zetaLP, IPC_zetaNotch, iStatus, IPC_ControlMode, Y_ControlMode, NumBl, IPC_PitComF) + IF ((IPC_ControlMode == 1) .OR. (Y_ControlMode == 2)) THEN + CALL IPC(rootMOOP, Azimuth, IPC_phi, Y_MErr, DT, IPC_KI, Y_IPC_KP, Y_IPC_KI, IPC_omegaHP, IPC_omegaLP, IPC_omegaNotch, IPC_zetaHP, IPC_zetaLP, IPC_zetaNotch, iStatus, IPC_ControlMode, Y_ControlMode, NumBl, IPC_PitComF) ELSE IPC_PitComF = 0.0 END IF @@ -575,7 +619,7 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM Y_ErrLPFFast = LPFilter(Y_MErr, DT, Y_omegaLPFast, iStatus, .FALSE., 2) ! Fast low pass filtered yaw error with a frequency of 1 Y_ErrLPFSlow = LPFilter(Y_MErr, DT, Y_omegaLPSlow, iStatus, .FALSE., 3) ! Slow low pass filtered yaw error with a frequency of 1/60 - Y_AccErr = Y_AccErr + ElapTime*SIGN(Y_ErrLPFFast**2, Y_ErrLPFFast) ! Integral of the fast low pass filtered yaw error + Y_AccErr = Y_AccErr + DT*SIGN(Y_ErrLPFFast**2, Y_ErrLPFFast) ! Integral of the fast low pass filtered yaw error IF (ABS(Y_AccErr) >= Y_ErrThresh) THEN ! Check if accumulated error surpasses the threshold Y_YawEndT = ABS(Y_ErrLPFSlow/Y_YawRate) + Time ! Yaw to compensate for the slow low pass filtered error @@ -591,7 +635,7 @@ SUBROUTINE DISCON (avrSWAP, aviFAIL, accINFILE, avcOUTNAME, avcMSG) BIND (C, NAM ! Output debugging information if requested: IF (DbgOut) THEN - WRITE (UnDb,FmtDat) Time, PitComT, PC_SpdErr, PC_KP, PC_KI, Y_MErr, rootMOOP(1), VS_RtPwr, GenTrq, testValue, PC_KPtest, PC_KItest + WRITE (UnDb,FmtDat) Time, PitComT, PC_SpdErr, PC_KP, PC_KI, Y_MErr, rootMOOP(1), VS_RtPwr, GenTrq, REAL(VS_ControlModeNow), REAL(PC_ControlModeNow) WRITE (UnDb2,FmtDat) Time, avrSWAP(1:85) END IF diff --git a/Source/FunctionToolbox.f90 b/Source/FunctionToolbox.f90 index dc10117f..ecb78684 100644 --- a/Source/FunctionToolbox.f90 +++ b/Source/FunctionToolbox.f90 @@ -40,7 +40,7 @@ REAL FUNCTION ratelimit(refSignal, measSignal, minRate, maxRate, DT) END FUNCTION ratelimit !------------------------------------------------------------------------------------------------------------------------------- ! PI controller, with output saturation - REAL FUNCTION PIController(error, kp, ki, minValue, maxValue, DT, I0, inst) + REAL FUNCTION PIController(error, kp, ki, minValue, maxValue, DT, I0, reset, inst) ! IMPLICIT NONE @@ -53,34 +53,33 @@ REAL FUNCTION PIController(error, kp, ki, minValue, maxValue, DT, I0, inst) REAL(4), INTENT(IN) :: DT INTEGER(4), INTENT(IN) :: inst REAL(4), INTENT(IN) :: I0 + LOGICAL, INTENT(IN) :: reset ! Local - INTEGER(4) :: i ! Counter for making arrays - REAL(4) :: PTerm ! Proportional term - REAL(4), DIMENSION(99), SAVE :: ITerm ! Integral term, current. - REAL(4), DIMENSION(99), SAVE :: ITermLast ! Integral term, the last time this controller was called. Supports 99 separate instances. - INTEGER(4), DIMENSION(99), SAVE :: FirstCall = (/ (1, i=1,99) /) ! First call of this function? + INTEGER(4) :: i ! Counter for making arrays + REAL(4) :: PTerm ! Proportional term + REAL(4), DIMENSION(99), SAVE :: ITerm = (/ (real(9999.9), i = 1,99) /) ! Integral term, current. + REAL(4), DIMENSION(99), SAVE :: ITermLast = (/ (real(9999.9), i = 1,99) /) ! Integral term, the last time this controller was called. Supports 99 separate instances. + INTEGER(4), DIMENSION(99), SAVE :: FirstCall = (/ (1, i=1,99) /) ! First call of this function? ! Initialize persistent variables/arrays, and set inital condition for integrator term - IF (FirstCall(inst) == 1) THEN - ITerm(1:99) = (/ (real(9999.9), i = 1,99) /) - ITermLast(1:99) = (/ (real(9999.9), i = 1,99) /) - + IF ((FirstCall(inst) == 1) .OR. reset) THEN ITerm(inst) = I0 ITermLast(inst) = I0 FirstCall(inst) = 0 - END IF + PIController = I0 + ELSE - PTerm = kp*error - ITerm(inst) = ITerm(inst) + DT*ki*error - ITerm(inst) = saturate(ITerm(inst), minValue, maxValue) - PIController = PTerm + ITerm(inst) - PIController = saturate(PIController, minValue, maxValue) - - ITermLast(inst) = ITerm(inst) + PTerm = kp*error + ITerm(inst) = ITerm(inst) + DT*ki*error + ITerm(inst) = saturate(ITerm(inst), minValue, maxValue) + PIController = PTerm + ITerm(inst) + PIController = saturate(PIController, minValue, maxValue) + + ITermLast(inst) = ITerm(inst) + END IF - ! WRITE (UnDb,FmtDat) real(FirstCall(inst)), ITerm(inst), error, kp, ki, DT END FUNCTION PIController !------------------------------------------------------------------------------------------------------------------------------- ! interp1 1-D interpolation (table lookup), xData and yData should be monotonically increasing diff --git a/Source/IPC.f90 b/Source/IPC.f90 index bf3185b3..7fc93aca 100644 --- a/Source/IPC.f90 +++ b/Source/IPC.f90 @@ -1,6 +1,6 @@ !------------------------------------------------------------------------------------------------------------------------------- ! Individual pitch control subroutine -SUBROUTINE IPC(rootMOOP, aziAngle, phi, Y_MErr, DT, KInter, omegaHP, omegaLP, omegaNotch, zetaHP, zetaLP, zetaNotch, iStatus, IPC_ControlMode, Y_ControlMode, NumBl, PitComIPCF) +SUBROUTINE IPC(rootMOOP, aziAngle, phi, Y_MErr, DT, KInter, Y_IPC_KP, Y_IPC_KI, omegaHP, omegaLP, omegaNotch, zetaHP, zetaLP, zetaNotch, iStatus, IPC_ControlMode, Y_ControlMode, NumBl, PitComIPCF) !............................................................................................................................... USE :: FunctionToolbox @@ -24,14 +24,16 @@ SUBROUTINE IPC(rootMOOP, aziAngle, phi, Y_MErr, DT, KInter, omegaHP, omegaLP, om REAL(4), INTENT(IN) :: omegaNotch ! Notch filter frequency REAL(4), INTENT(IN) :: phi ! Phase offset added to the azimuth angle REAL(4), INTENT(IN) :: rootMOOP(3) ! Root out of plane bending moments of each blade - REAL(4), INTENT(IN) :: Y_MErr ! Yaw alignment error [rad] + REAL(4), INTENT(IN) :: Y_MErr ! Yaw alignment error, measured [rad] REAL(4), INTENT(IN) :: zetaHP ! High-pass filter damping value REAL(4), INTENT(IN) :: zetaLP ! Low-pass filter damping value REAL(4), INTENT(IN) :: zetaNotch ! Notch filter damping value INTEGER(4), 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(4), INTENT(IN) :: NumBl ! Number of turbine blades INTEGER(4), INTENT(IN) :: Y_ControlMode ! Yaw control mode - INTEGER(4), INTENT(IN) :: IPC_ControlMode ! Turn Individual Pitch Control (IPC) for fatigue load reductions (pitch contribution) on = 1/off = 1 + REAL(4), INTENT(IN) :: Y_IPC_KP ! Yaw-by-IPC proportional controller gain Kp + REAL(4), INTENT(IN) :: Y_IPC_KI ! Yaw-by-IPC integral controller gain Ki + INTEGER(4), INTENT(IN) :: IPC_ControlMode ! Turn Individual Pitch Control (IPC) for fatigue load reductions (pitch contribution) on = 1/off = 0 ! Outputs @@ -59,7 +61,7 @@ SUBROUTINE IPC(rootMOOP, aziAngle, phi, Y_MErr, DT, KInter, omegaHP, omegaLP, om ! Calculate commanded IPC pitch angles - CALL CalculatePitCom(rootMOOPF, aziAngle, Y_MErr, DT, KInter, omegaHP, zetaHP, phi, iStatus, IPC_ControlMode, Y_ControlMode, PitComIPC) + CALL CalculatePitCom(rootMOOPF, aziAngle, Y_MErr, DT, KInter, Y_IPC_KP, Y_IPC_KI, omegaHP, zetaHP, phi, iStatus, IPC_ControlMode, Y_ControlMode, PitComIPC) ! Filter PitComIPC with second order low pass filter @@ -80,7 +82,7 @@ SUBROUTINE IPC(rootMOOP, aziAngle, phi, Y_MErr, DT, KInter, omegaHP, omegaLP, om ! Calculates the commanded pitch angles. ! NOTE: if it is required for this subroutine to be used multiple times (for 1p and 2p IPC for example), the saved variables ! IntAxisTilt and IntAxisYaw need to be modified so that they support multiple instances (see LPFilter in the Filters module). - SUBROUTINE CalculatePitCom(rootMOOP, aziAngle, Y_MErr, DT, KInter, omegaHP, zetaHP, phi, iStatus, IPC_ControlMode, Y_ControlMode, PitComIPC) + SUBROUTINE CalculatePitCom(rootMOOP, aziAngle, Y_MErr, DT, KInter, Y_IPC_KP, Y_IPC_KI, omegaHP, zetaHP, phi, iStatus, IPC_ControlMode, Y_ControlMode, PitComIPC) !............................................................................................................................... IMPLICIT NONE @@ -96,7 +98,10 @@ SUBROUTINE CalculatePitCom(rootMOOP, aziAngle, Y_MErr, DT, KInter, omegaHP, zeta REAL(4), INTENT(IN) :: zetaHP ! High-pass filter damping value INTEGER(4), 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(4), INTENT(IN) :: Y_ControlMode ! Yaw control mode - INTEGER(4), INTENT(IN) :: IPC_ControlMode ! Turn Individual Pitch Control (IPC) for fatigue load reductions (pitch contribution) on = 1/off = 1 + REAL(4), INTENT(IN) :: Y_MErr ! Yaw alignment error, measured [rad] + REAL(4), INTENT(IN) :: Y_IPC_KP ! Yaw-by-IPC proportional controller gain Kp + REAL(4), INTENT(IN) :: Y_IPC_KI ! Yaw-by-IPC integral controller gain Ki + INTEGER(4), INTENT(IN) :: IPC_ControlMode ! Turn Individual Pitch Control (IPC) for fatigue load reductions (pitch contribution) on = 1/off = 0 ! Outputs @@ -107,7 +112,7 @@ SUBROUTINE CalculatePitCom(rootMOOP, aziAngle, Y_MErr, DT, KInter, omegaHP, zeta REAL(4) :: axisTilt, axisYaw, axisYawF ! Direct axis and quadrature axis outputted by Coleman transform REAL(4), SAVE :: IntAxisTilt, IntAxisYaw ! Integral of the direct axis and quadrature axis REAL(4) :: IntAxisYawIPC ! IPC contribution with yaw-by-IPC component - REAL(4) :: Y_MErr, Y_MErrF, Y_MErrF_IPC ! Unfiltered and filtered yaw alignment error [rad] + REAL(4) :: Y_MErrF, Y_MErrF_IPC ! Unfiltered and filtered yaw alignment error [rad] REAL(4) :: PitComIPC_woYaw(3) ! Debugging @@ -144,7 +149,7 @@ SUBROUTINE CalculatePitCom(rootMOOP, aziAngle, Y_MErr, DT, KInter, omegaHP, zeta IF (Y_ControlMode == 2) THEN axisYawF = HPFilter(axisYaw, DT, omegaHP, iStatus, 1) Y_MErrF = SecLPFilter(Y_MErr, DT, omegaLP, zetaLP, iStatus, 1) - Y_MErrF_IPC = PIController(Y_MErrF, -0.16, -0.002, -100.0, 100.0, DT, 0.0, 3) + Y_MErrF_IPC = PIController(Y_MErrF, Y_IPC_KP, Y_IPC_KI, -100.0, 100.0, DT, 0.0, .FALSE., 3) ELSE axisYawF = axisYaw Y_MErrF = 0.0