Unverified Commit cadf5cbd authored by David Verelst's avatar David Verelst Committed by GitHub
Browse files

Merge pull request #29 from DTUWindEnergy/check_iteration

Check iteration
parents 3daf7ffd 3a9d90cf
......@@ -10,6 +10,7 @@ module dtu_we_controller
integer CtrlStatus
real(mk) dump_array(50)
real(mk) time_old
logical repeated
contains
!**************************************************************************************************
subroutine init_regulation(array1, array2) bind(c, name='init_regulation')
......@@ -278,6 +279,7 @@ subroutine init_regulation(array1, array2) bind(c, name='init_regulation')
! Initiate the dynamic variables
stepno = 0
time_old = 0.0_mk
repeated=.FALSE.
AddedPitchRate = 0.0_mk
PitchAngles=0.0_mk
AveragedMeanPitchAngles=0.0_mk
......@@ -429,7 +431,6 @@ subroutine update_regulation(array1, array2) bind(c,name='update_regulation')
! Local variables
integer GridFlag, EmergPitchStop, ActiveMechBrake
real(mk) GenSpeed, wsp, PitchVect(3), Pe, TT_acc(2), time
real(mk) GenTorqueRef, PitchColRef
EmergPitchStop = 0
ActiveMechBrake = 0
! Time
......@@ -437,10 +438,20 @@ subroutine update_regulation(array1, array2) bind(c,name='update_regulation')
!***********************************************************************************************
! Increment time step (may actually not be necessary in type2 DLLs)
!***********************************************************************************************
!somehow the controller gets called twice in the very first time step...
if ((time==deltat).AND. (repeated==.FALSE.)) then
time_old=0.0_mk
repeated=.TRUE.
endif
if (time .gt. time_old) then
deltat = time - time_old
time_old = time
stepno = stepno + 1
newtimestep = .TRUE.
PitchColRefOld = PitchColRef
GenTorqueRefOld = GenTorqueRef
else
newtimestep = .FALSE.
endif
! Rotor (Generator) speed in LSS
GenSpeed = array1(2)/GearRatio
......
......@@ -8,15 +8,16 @@ module dtu_we_controller_fcns
! Constants
integer maxwplines
parameter(maxwplines=100)
logical :: newtimestep=.TRUE.
! Types
type Tpidvar
real(mk) Kpro,Kdif,Kint,outmin,outmax,velmax,error1,outset1,outres1
integer stepno1
integer :: stepno1 = 0
real(mk) outset,outpro,outdif,error1_old,outset1_old,outres1_old,outres
end type Tpidvar
type Tpid2var
real(mk) Kpro(2),Kdif(2),Kint(2),outmin,outmax,velmax,error1(2),outset1,outres1
integer stepno1
integer :: stepno1 = 0
real(mk) outset,outpro,outdif,error1_old(2),outset1_old,outres1_old,outres
end type Tpid2var
type Twpdata
......@@ -140,9 +141,9 @@ function PID(stepno, dt, kgain, PIDvar, error)
PIDvar%error1_old = PIDvar%error1
endif
! Update the integral term
PIDvar%outset = PIDvar%outset1_old + 0.5_mk*(error + PIDvar%error1)*Kgain(2)*PIDvar%Kint*dt
PIDvar%outset = PIDvar%outset1_old + 0.5_mk*(error + PIDvar%error1_old)*Kgain(2)*PIDvar%Kint*dt
! Update proportional term
PIDvar%outpro = Kgain(1)*PIDvar%Kpro*0.5_mk*(error + PIDvar%error1)
PIDvar%outpro = Kgain(1)*PIDvar%Kpro*0.5_mk*(error + PIDvar%error1_old)
! Update differential term
PIDvar%outdif = Kgain(3)*PIDvar%Kdif*(error - PIDvar%error1_old)/dt
! Sum to up
......@@ -198,11 +199,11 @@ function PID2(stepno,dt,kgain,PIDvar,error,added_term)
PIDvar%error1_old = PIDvar%error1
endif
! Update the integral term
PIDvar%outset = PIDvar%outset1_old + 0.5_mk*dt*(Kgain(2, 1)*PIDvar%Kint(1)*(error(1) + PIDvar%error1(1))&
+Kgain(2, 2)*PIDvar%Kint(2)*(error(2) + PIDvar%error1(2)))
PIDvar%outset = PIDvar%outset1_old + 0.5_mk*dt*(Kgain(2, 1)*PIDvar%Kint(1)*(error(1) + PIDvar%error1_old(1))&
+Kgain(2, 2)*PIDvar%Kint(2)*(error(2) + PIDvar%error1_old(2)))
! Update proportional term
PIDvar%outpro = 0.5_mk*(Kgain(1, 1)*PIDvar%Kpro(1)*(error(1) + PIDvar%error1(1))&
+Kgain(1, 2)*PIDvar%Kpro(2)*(error(2) + PIDvar%error1(2)))
PIDvar%outpro = 0.5_mk*(Kgain(1, 1)*PIDvar%Kpro(1)*(error(1) + PIDvar%error1_old(1))&
+Kgain(1, 2)*PIDvar%Kpro(2)*(error(2) + PIDvar%error1_old(2)))
! Update differential term
PIDvar%outdif = (Kgain(3, 1)*PIDvar%Kdif(1)*(error(1) - PIDvar%error1_old(1)))/dt + added_term*dt
! Sum up
......
......@@ -11,31 +11,32 @@ module misc_mod
! First order filter
type Tfirstordervar
real(mk) tau, x1, x1_old, y1, y1_old
integer stepno1
integer :: stepno1 = 0
end type Tfirstordervar
! Second order low pass filter filter
type Tlowpass2order
real(mk) zeta, f0, x1, x2, x1_old, x2_old, y1, y2, y1_old, y2_old
integer stepno1
integer :: stepno1 = 0
end type Tlowpass2order
! Second order notch filter
type Tnotch2order
real(mk) :: zeta1 = 0.1_mk
real(mk) :: zeta2 = 0.001_mk
real(mk) f0, x1, x2, x1_old, x2_old, y1, y2, y1_old, y2_old
integer stepno1
integer :: stepno1 = 0
end type Tnotch2order
! Second order band-pass filter
type Tbandpassfilt
real(mk) :: zeta = 0.02_mk
real(mk) :: tau = 0.0_mk
real(mk) f0, x1, x2, x1_old, x2_old, y1, y2, y1_old, y2_old
integer stepno1
integer :: stepno1 = 0
end type Tbandpassfilt
! Time delay
type Ttdelay
real(mk) xz(40)
integer stepno1
real(mk) xz_old(40)
integer :: stepno1 = 0
end type Ttdelay
contains
!**************************************************************************************************
......@@ -105,9 +106,9 @@ function lowpass2orderfilt(dt, stepno, filt, x)
y = a1*filt%y1_old + a2*filt%y2_old + b0*x + b1*filt%x1_old + b2*filt%x2_old
endif
! Save previous values
filt%x2 = filt%x1
filt%x2 = filt%x1_old
filt%x1 = x
filt%y2 = filt%y1
filt%y2 = filt%y1_old
filt%y1 = y
filt%stepno1=stepno
! Output
......@@ -153,9 +154,9 @@ function notch2orderfilt(dt,stepno,filt,x)
y = a1*filt%y1_old + a2*filt%y2_old + b0*x + b1*filt%x1_old + b2*filt%x2_old
endif
! Save previous values
filt%x2 = filt%x1
filt%x2 = filt%x1_old
filt%x1 = x
filt%y2 = filt%y1
filt%y2 = filt%y1_old
filt%y1 = y
filt%stepno1 = stepno
! Output
......@@ -200,9 +201,9 @@ function bandpassfilt(dt, stepno, filt, x)
y = a1*filt%y1_old + a2*filt%y2_old + b0*x + b1*filt%x1_old + b2*filt%x2_old
endif
! Save previous values
filt%x2 = filt%x1
filt%x2 = filt%x1_old
filt%x1 = x
filt%y2 = filt%y1
filt%y2 = filt%y1_old
filt%y1 = y
filt%stepno1 = stepno
! Output
......@@ -222,8 +223,11 @@ function timedelay(dt, stepno, filt, Td, x)
filt%xz(k) = x
end do
endif
if (stepno .gt. filt%stepno1) then
filt%xz_old = filt%xz
endif
do k = 40, 2, -1
filt%xz(k) = filt%xz(k - 1)
filt%xz(k) = filt%xz_old(k - 1)
end do
filt%xz(1) = x
! Output
......@@ -232,6 +236,7 @@ function timedelay(dt, stepno, filt, Td, x)
else
timedelay = filt%xz(n)
endif
filt%stepno1 = stepno
return
end function timedelay
!**************************************************************************************************
......
......@@ -18,8 +18,16 @@ module turbine_controller_mod
real(mk) TAve_Pitch,DeltaPitchThreshold,AveragedMeanPitchAngles(3),AveragedPitchReference(3)
! Dynamic variables
integer :: stepno = 0, w_region = 0
real(mk) AddedPitchRate, PitchColRef0, GenTorqueRef0, PitchColRefOld, GenTorqueRefOld
real(mk) TimerGenCutin, TimerStartup, TimerExcl, TimerShutdown, TimerShutdown2
real(mk) AddedPitchRate, PitchColRef0, GenTorqueRef0
real(mk) :: PitchColRef = 0.0_mk
real(mk) :: GenTorqueRef = 0.0_mk
real(mk) :: PitchColRefOld = 0.0_mk
real(mk) :: GenTorqueRefOld = 0.0_mk
real(mk) :: TimerGenCutin = 0.0_mk
real(mk) :: TimerStartup = 0.0_mk
real(mk) :: TimerExcl = 0.0_mk
real(mk) :: TimerShutdown = 0.0_mk
real(mk) :: TimerShutdown2 = 0.0_mk
real(mk) GenSpeed_at_stop, GenTorque_at_stop
real(mk) excl_flag
real(mk)::outmax_old=0.0_mk,outmin_old=0.0_mk
......@@ -94,11 +102,6 @@ subroutine turbine_controller(CtrlStatus, GridFlag, GenSpeed, PitchVect, wsp, Pe
if (CtrlStatus .gt. 0) then
call shut_down(CtrlStatus, GenSpeed, PitchVect, wsp, GenTorqueRef, PitchColRef, dump_array)
endif
!***********************************************************************************************
! Save reference signals
!***********************************************************************************************
PitchColRefOld = PitchColRef
GenTorqueRefOld = GenTorqueRef
return
end subroutine turbine_controller
!**************************************************************************************************
......@@ -158,7 +161,7 @@ subroutine normal_operation(GenSpeed, PitchVect, wsp, Pe, TTfa_acc, GenTorqueRef
! Active DT damping based on filtered rotor speed
!***********************************************************************************************
call drivetraindamper(GenSpeed, Qdamp_ref, dump_array)
TimerGenCutin = TimerGenCutin + deltat
if (newtimestep) TimerGenCutin = TimerGenCutin + deltat
x = switch_spline(TimerGenCutin, CutinVar%delay, 2.0_mk*CutinVar%delay)
GenTorqueRef = min(max(GenTorqueRef + Qdamp_ref*x, 0.0_mk), GenTorqueMax)
!***********************************************************************************************
......@@ -244,7 +247,7 @@ subroutine start_up(CtrlStatus, GenSpeed, PitchVect, wsp, GenTorqueRef, PitchCol
dummy = PID(stepno, deltat, kgain_torque, PID_gen_var, GenSpeedFiltErr)
elseif (TimerStartup.lt.CutinVar%delay) then
! Start increasing the timer for the delay
TimerStartup = TimerStartup + deltat
if (newtimestep) TimerStartup = TimerStartup + deltat
! Generator is cut-in
generator_cutin = .true.
! Gradually set the minimum pitch angle to optimal pitch
......@@ -309,14 +312,14 @@ subroutine shut_down(CtrlStatus, GenSpeed, PitchVect, wsp, GenTorqueRef, PitchCo
PitchMeanFilt = lowpass1orderfilt(deltat, stepno, pitchfirstordervar, PitchMean)
PitchMeanFilt = min(PitchMeanFilt, 30.0_mk*degrad) ! Maximum of 30 deg
! Start increasing the timer for the delay
TimerShutdown2 = TimerShutdown2 + deltat
if (newtimestep) TimerShutdown2 = TimerShutdown2 + deltat
! Generator settings
select case(CtrlStatus)
case(1, 3, 4, 7)
if (GenSpeed .gt. GenSpeed_at_stop*0.8_mk) then
GenTorqueRef = GenTorque_at_stop
else
TimerShutdown = TimerShutdown + deltat
if (newtimestep) TimerShutdown = TimerShutdown + deltat
GenTorqueRef = max(0.0_mk, GenTorque_at_stop*(1.0_mk - TimerShutdown/CutoutVar%torquedelay))
endif
case(2, 5, 6)
......@@ -339,7 +342,7 @@ subroutine shut_down(CtrlStatus, GenSpeed, PitchVect, wsp, GenTorqueRef, PitchCo
deltat*max(CutoutVar%pitchvelmax, CutoutVar%pitchvelmax2))
case(-2) ! Pitch-out before cut-in
PitchColRef = min(PitchStopAng, PitchColRefOld + deltat*CutoutVar%pitchvelmax)
end select
end select
! Write into dump array
dump_array(1) = GenTorqueRef*GenSpeed
dump_array(3) = y(1)
......@@ -349,7 +352,7 @@ subroutine shut_down(CtrlStatus, GenSpeed, PitchVect, wsp, GenTorqueRef, PitchCo
return
end subroutine shut_down
!**************************************************************************************************
subroutine monitoring(CtrlStatus, GridFlag, GenSpeed, TTAcc, PitchVect, PitchColRef, dump_array)
subroutine monitoring(CtrlStatus, GridFlag, GenSpeed, TTAcc, PitchVect, PitchColRefmonitor, dump_array)
!
! Lower level system monitoring. It changes the controller status to:
! - (1) if filtered GenSpeed is higher than the overspeed limit.
......@@ -364,7 +367,7 @@ subroutine monitoring(CtrlStatus, GridFlag, GenSpeed, TTAcc, PitchVect, PitchCol
real(mk), intent(in) :: TTAcc ! Tower top acceleration [m/s**2].
real(mk), intent(in) :: GenSpeed ! Measured generator speed [rad/s].
real(mk), intent(in) :: PitchVect(3) ! Measured pitch angles [rad].
real(mk), intent(in) :: PitchColRef ! Reference collective pitch [rad].
real(mk), intent(in) :: PitchColRefmonitor ! Reference collective pitch [rad].
real(mk), intent(inout) :: dump_array(50) ! Array for output.
real(mk) GenSpeedFilt, dGenSpeed_dtFilt, TTAccFilt
real(mk) y(2),DiffPitch
......@@ -429,9 +432,9 @@ subroutine monitoring(CtrlStatus, GridFlag, GenSpeed, TTAcc, PitchVect, PitchCol
PitchRefs(2:NAve_Pitch,1:3)=PitchRefs(1:NAve_Pitch-1,1:3)
! Update new values
PitchAngles(1,1:3)=PitchVect
PitchRefs(1,1)=PitchColRef
PitchRefs(1,2)=PitchColRef
PitchRefs(1,3)=PitchColRef
PitchRefs(1,1)=PitchColRefmonitor
PitchRefs(1,2)=PitchColRefmonitor
PitchRefs(1,3)=PitchColRefmonitor
! Add new values
AveragedMeanPitchAngles=AveragedMeanPitchAngles+PitchAngles(1,1:3)/dfloat(NAve_Pitch)
AveragedPitchReference=AveragedPitchReference+PitchRefs(1,1:3)/dfloat(NAve_Pitch)
......@@ -479,6 +482,10 @@ subroutine torquecontroller(GenSpeed, GenSpeedFilt, dGenSpeed_dtFilt, PitchMean,
!***********************************************************************************************
! Speed ref. changes max. <-> min. for torque contr. and remains at rated for pitch contr.
!***********************************************************************************************
if (newtimestep) then
outmax_old=PID_gen_var%outmax
outmin_old=PID_gen_var%outmin
endif
select case (PartialLoadControlMode)
case (1)
if (GenSpeedFilt .gt. 0.5_mk*(GenSpeedRefMax + GenSpeedRefMin)) then
......@@ -552,8 +559,6 @@ subroutine torquecontroller(GenSpeed, GenSpeedFilt, dGenSpeed_dtFilt, PitchMean,
if ((abs(outmin-outmin_old)/deltat) .gt. PID_gen_var%velmax) then
outmin = outmin_old + dsign(PID_gen_var%velmax*deltat, outmin-outmin_old)
endif
outmax_old=outmax
outmin_old=outmin
PID_gen_var%outmin = outmin
PID_gen_var%outmax = outmax
if (PID_gen_var%outmin .gt. PID_gen_var%outmax) PID_gen_var%outmin = PID_gen_var%outmax
......@@ -660,7 +665,7 @@ subroutine rotorspeedexcl(GenSpeedFilt, GenTorque, Qg_min_partial, GenTorqueMax_
x = 0.0_mk
y = notch2orderfilt(deltat, stepno, ExcluZone%notch, GenSpeedFilt)
GenSpeedFiltNotch = y(1)
TimerExcl = TimerExcl + deltat
if (newtimestep) TimerExcl = TimerExcl + deltat
select case (w_region)
case (0)
if ((GenSpeedFilt .gt. Lwr*0.99_mk) .and. (TimerGenCutin .gt. CutinVar%delay))then
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment