R/ZTTrainTheoTurbineModel.R

In AndyClifton/PowerPerformance: Package to analyze wind resource and wind turbine data (developer test version)

#' Initialize the zero-turbulence power curve
#' 
#' \code{ZTTrainTheoTurbineModel} calculates the theoretical zero-turbulence power curve
#' using the method set out in the draft version of IEC 61400-12-1 (2013). This 
#' should be the second stage in calculating a site-specific power curve using 
#' the zero-turbulence power curve methdology.
#' 
#' The input to this function is generated by \code{\link{ZTTrainInitTurbineModel}}.
#' 
#' The result of this function is used as the input to \code{\link{ZTTrainFinalTurbineModel}}.
#' 
#' @param PC.param a data.frame of parameters describing the initial
#'   zero-turbulence power curve, created using \code{\link{ZTTrainInitTurbineModel}}
#' @param ws bin-averaged wind speeds
#' @param Ti bin-averaged turbulence intensity
#' @param power bin-averaged power
#' @param rho density
#' @return PC.values a data.frame containing binned wind speed and power values
#' @export
#' 
#' @family zero-turbulence power curve methods
ZTTrainTheoTurbineModel <- function(PC.param,
                                    ws,
                                    Ti,
                                    power,
                                    rho = 1.225,
                                    limit.p.sim = TRUE){
  
  # First get the theoretical power curve
  PC.values = data.frame('ws.binmean' = ws,
                         'ti.binmean' = Ti)
  
  ## START OF ITERATION
  # set defaults for testing
  test.converged = FALSE
  
  while (test.converged == FALSE){
    # iterative loop
    # note that we want the most recent estimate of the power curve parameters
    iter = nrow(PC.param)
    # get the theoretical power curve. 
    ###################################################
    ## This appears to be where Rozenn and I differ. ##
    ###################################################
    PC.values$power.binmean <- PCfromParam(data.frame(rho = rho,
                                                      ws = PC.values$ws.binmean,
                                                      cp = PC.param$cp.max[iter],
                                                      ws.cutin = PC.param$ws.cutin[iter],
                                                      ws.rated = PC.param$ws.rated[iter],
                                                      power.rated = PC.param$power.rated[iter],
                                                      diameter = PC.param$diameter[1]))[,"power.binmean"]
    # remove empty bins
    PC.values <- PC.values[!is.na(PC.values$ti.binmean),]    
    
    # Get the simulated power (kW) in 10 minutes using the turbulence intensity ----    
    sim.power <- SimPC(PC.param = PC.param[iter,],
                       ws = PC.values$ws.binmean,
                       Ti = PC.values$ti.binmean,
                       rho = rho)
        
    # now get the new values for the power curve
    sim.P.max = max(sim.power, na.rm = TRUE)
    sim.ws.cutin = PC.values$ws.binmean[min(which(sim.power >= (0.1/100)*PC.param$power.rated[iter]))]
    # get cp
    sim.cp = 2*sim.power / ((rho * PC.values$ws.binmean^3 * (pi * PC.param$diameter[1]^2/4))/1000)
    
    sim.cp.max = max(sim.cp, na.rm = TRUE)
    sim.ws.rated = (2*sim.P.max*1000 / (rho * sim.cp.max * pi * PC.param$diameter^2/4))^(1/3)
    
    # GET VALUES FOR NEXT ITERATION
    # default assumption is they don't change.
    new.power.rated = PC.param$power.rated[iter]
    new.ws.cutin = PC.param$ws.cutin[iter]
    new.cp.max = PC.param$cp.max[iter]
    new.ws.rated = PC.param$ws.rated[iter]
    
    # CHECK RESULTS
    # assume we passed, check to see if that's true
    # depending on what converged (or not), change certain things    
    if ((abs(sim.P.max - PC.param$power.rated[1]) >= ((0.1/100)*PC.param$power.rated[1]))){
      # get the new value of power.rated      
      new.power.rated = PC.param$power.rated[iter] - sim.P.max + PC.param$power.rated[1]
    } else {
      if (abs(sim.ws.cutin - PC.param$ws.cutin[1]) >= 0.5){
        # get the new value of ws.cutin
        new.ws.cutin = PC.param$ws.cutin[iter] - sim.ws.cutin + PC.param$ws.cutin[1]
      } else {        
        # power and cut in are correct
        if (abs(sim.cp.max - PC.param$cp.max[1]) >= 0.01){
          # get the new value of cp.max
          new.cp.max = PC.param$cp.max[iter] - sim.cp.max + PC.param$cp.max[1]          
        } else {
          test.converged = TRUE
        }
      }
    }
    # deadman's switch
    if (iter > 20){
      test.converged = TRUE
    }
    
    new.ws.rated = ((2 * new.power.rated*1000) /
                      (rho * new.cp.max * pi * PC.param$diameter[1]^2/4))^(1/3)
    
    # INITIAL VALUES FOR NEXT ITERATION
    newvals <- data.frame("iteration" = iter,
                          "power.rated" = new.power.rated,
                          "ws.cutin" = new.ws.cutin,                          
                          "cp.max" = new.cp.max,
                          "ws.rated" = new.ws.rated,
                          "diameter" = PC.param$diameter[1])    
    PC.param <- rbind(PC.param,
                      newvals)
    
    # THEORETICAL POWER CURVES
    
  }
  # END OF ITERATION
  # return the parameters that define the theoretical zero turbulence power cuver
  return(list("values" = PC.values,
              "param" = newvals,
              "iterations" = PC.param))
}