R [Shiny] How to implement helper functions into a dynamical model using deSolve and Shiny

Question

I am trying to implement a Vensim Model that I built, using R-Shiny and the deSolve Package. The model itself is inspired by a Report from Rethink X but I tried to adapt it to the swiss situation.

It basically simulates an energy grid which is solely powered by Wind, Solar, Li-Batteries, and Hydro Power.

I've had trouble getting the model running in R-Shiny as I have never used a helper function in combination with a dynamic model in R.

Here is the code:

library(shiny)
library(deSolve)
library(ggplot2)
library(gridExtra)

#TODO define vectors for the data inputs
set.seed(1)
solarData <- rnorm(FINISH, mean = 0.5,sd=1)
set.seed(1)
windData <- rnorm(FINISH, mean = 0.5,sd=1)
set.seed(1)
demandData <- rnorm(FINISH, mean = 0.5,sd=1)
set.seed(1)
hydroRunOffData <- rnorm(FINISH, mean = 0.5,sd=1)


ui <- fluidPage(
  sliderInput("iSolarCap", "Solar Capacity", min = 20000, max = 150000, step = 2000, value = 50000),
  sliderInput("iWindCap", "Wind Capacity", min = 50, max = 500, step = 50, value = 50),
  sliderInput("iGrowthFactorDemand", "Growth Factor Demand", min = 1, max = 1.5, step = 0.05, value = 1),
  sliderInput("iLithiumCap", "Lithium Capacity", min = 1000, max = 100000, step = 1000, value = 10000),
  
  plotOutput(outputId = "arrange")
)



server <- function(input, output, session) {
  
  #defining Values that can be changed
  solarCap <- reactiveVal(1)
  windCap <- reactiveVal(1)
  growthFactorDemand <- reactiveVal(1)
  lithiumCap <- reactiveVal(1)

  #creating simulation timeline vector
  START <-1; FINISH<-17522; STEP<-1
  simtime <- seq(START, FINISH, by = STEP)
  
  #TODO define vectors for the data inputs
  solarData <- soldata
  windData <- windata
  demandData <- demdata
  hydroRunOffData <- rivdata


  
  observe({
    
    #linking input with functions ??not sure if thats what it does
    solarCap(input$iSolarCap)
    windCap(input$iWindCap)
    growthFactorDemand(input$iGrowthFactorDemand)
    lithiumCap(input$iLithiumCap)
    
    #Assigning Stocks
    stocks   <- c(sLithium = lithiumCap(),
                  sPumpedHydro =1,
                  sSuperPower=0)
    

    
    #Assigning Auxiliary variables that are constants
    auxs <- list(aSolarCap = solarCap(),
                 aWindCap  = windCap(),
                 aGrowthFactorDemand = growthFactorDemand(),
                 aEfficiencyLithium = 0.95,
                 aCLithium = 0.5,
                 aLithiumCap = lithiumCap(),
                 aEfficiencyPumpedHydro = 0.75,
                 aChargePowerPumHydro = 5000,
                 aPumpedHydroCap = 150000)
    
    #Maybe here is the spot to put the various if Then Functions ????
    chargeFunction <- function(supplydemandGap){
      if (supplydemandGap>=0){
        return(supplydemandGap)
      }
    }
    
    dischargeFunction <- function(supplydemandGap){
      if (supplydemandGap<0){
        return(supplydemandGap*-1)
      }
    }
    
    #This Function checks how much Charge Lithium should get
    currentLithiumChargePower <- function(sLithium, aCharge, aLithiumCap, aEfficiencyLithium) {
      availableCap <- (aLithiumCap-sLithium)
      potCharge    <- (aCharge/aEfficiencyLithium)
      chargePower  <- (aCLithium*aLithiumCap)/aEfficiencyLithium
      
      if (sLithium < aLithiumCap){#check if there is still capacity left
        if (availableCap<potCharge){#check if available capacity is smaller than potential charge
          if (chargePower<=potCharge){#check if charge power is smaller than potential charge
            return(chargePower) #returning charge power if pot power higher
          }
          else{#returning potential charge if charge power not achieved
            return(potCharge)
          }
        }
        else{#if potential charge is higher than available capacity
          if (chargePower<=availableCap){#check if charge power is smaller than available Capacity
            return(chargePower) #returning charge power if available Capacity higher
          }
          else{
            return(availableCap)
          }
        }
      }
      else{
        print("Lithium is Full")
      }
    }
    
    #Checks how much charge pumped hydro should get
    currentPumHydroChargePower<- function(sPumpedHydro,aCharge,aLithiumCharging,aPumpedHydroCap,aChargePowerPumHydro){
      availableCap <- (aPumpedHydroCap-sPumpedHydro)
      potCharge    <- (aCharge-aLithiumCharging)/aEfficiencyPumpedHydro
      chargePower  <- aChargePowerPumHydro/aEfficiencyPumpedHydro
      
      if (sLithium < aLithiumCap){#check if there is still capacity left
        if (availableCap<potCharge){#check if available capacity is smaller than potential charge
          if (chargePower<=potCharge){#check if charge power is smaller than potential charge
            return(chargePower) #returning charge power if pot power higher
          }
          else{#returning potential charge if charge power not achieved
            return(potCharge)
          }
        }
        else{#if potential charge is higher than available capacity
          if (chargePower<=availableCap){#check if charge power is smaller than available Capacity
            return(chargePower) #returning charge power if available Capacity higher
          }
          else{
            return(availableCap)
          }
        }
      }
      else{
        print("Pumped Hydro is Full")
      }
    }
    
    #checks how much excess power is generated 
    currentExcessPower <- function(aCharge, aLithiumCharging, aPumHydroCharging){
      excessPower <- aCharge - aLithiumCharging - aPumpHydroCharging
      if (excessPower > 0){
        return(excessPower)
      }
    }
    
    #checks how much lithium can be discharged
    curLithiumDischarge <- function(sLithium, aDischarge,aCLithium,aLithiumCap){
      availableCap <- (sLithium)
      potDischarge <- (aDischarge)
      dischargePow <- (aCLithium*aLithiumCap)
      
      if (sLithium>0){#check if there is still capacity left
        if (availableCap>potDischarge){#check if available capacity is bigger than potential discharge
          if (dischargePow<potDischarge){#check if discharge power is smaller than potential discharge
            return(dischargePow) #returning charge power if pot power higher
          }
          else{#returning potential discharge if discharge power not achieved
            return(potDischarge)
          }
        }
        else{#if potential discharge is higher than available capacity (almost empty)
          if (dischargePow<=availableCap){#check if discharge power is smaller than available Capacity
            return(dischargePow) #returning charge power if available Capacity higher
          }
          else{
            return(availableCap)
          }
        }
      }
      else{
        print("Lithium is Empty")
      }
      
    }
    
    #checks how much Pumped Hydro can be discharged
    curPumHyDischarge <- function(sPumpedHydro, aDischarge,aLithiumDischarge,aChargePowerPumHydro){
      availableCap <- (sPumpedHydro)
      potDischarge <- (aDischarge- aLithiumDischarge)
      dischargePow <- (aChargePowerPumHydro)
      
      if (sPumpedHydro>0){#check if there is still capacity left
        if (availableCap>potDischarge){#check if available capacity is bigger than potential discharge
          if (dischargePow<potDischarge){#check if discharge power is smaller than potential discharge
            return(dischargePow) #returning charge power if pot power higher
          }
          else{#returning potential discharge if discharge power not achieved
            return(potDischarge)
          }
        }
        else{#if potential discharge is higher than available capacity (almost empty)
          if (dischargePow<=availableCap){#check if discharge power is smaller than available Capacity
            return(dischargePow) #returning charge power if available Capacity higher
          }
          else{
            return(availableCap)
          }
        }
      }
      else{
        print("Pumped Hydro is Empty")
      }
    }
    
    model <- function(time,stocks,auxs){
      with(as.list(c(stocks,auxs)),{
        #Getting the current input values
        aSolarProduction <- solarData[time]*aSolarCap
        aWindProduction  <- windData[time]*aWindCap
        aHydroRunOff     <- hydroRunOffData[time]
        aDemand          <- demandData[time]*aGrowthFactorDemand

        
        
        aSupplyDemandGap <- aDemand - (aSolarProduction+aWindProduction+aHydroRunOff)
        
        #checking if grid needs charge or discharge
        aCharge          <- chargeFunction(aSupplyDemandGap)
        aDischarge       <- dischargeFunction(aSupplyDemandGap)
        
        #Finding the amount that lithium can be charged
        aLithiumCharging <- currentLithiumChargePower(sLithium,aCharge, aLithiumCap,aEfficiencyLithium)
        fChargingLithium <- aLithiumCharging*aEfficiencyLithium
        
        #Finding the amount that Pumped Hydro can be charged
        aPumHydroCharging<- currentPumHydroChargePower(sPumpedHydro,aCharge,aLithiumCharging,aPumpedHydroCap,aChargePowerPumHydro)
        fChargingPumHydro<- aPumHydroCharging*aEfficiencyPumpedHydro
        
        #Excess Power Deliverance
        fSPCharging      <- currentExcessPower(aCharge, aLithiumCharging, aPumHydroCharging)
        
        #Finding discharge of Lithium
        aLithiumDischarge<- curLithiumDischarge(sLithium, aDischarge,aCLithium,aLithiumCap)
        fDischargingLi   <- aLithiumDischarge
        
        #Finding discharge of Pumped Hydro
        aPumHydroDisch   <- curPumHyDischarge(sPumpedHydro, aDischarge, aLithiumDischarge,aChargePowerPumHydro)
        fDischargingPumHy<- aPumHydroDisch
        
        
        
        ###Stored Hydro to be implemented

        
        dL_dt            <- fChargingLithium - fDischargingLi
        dP_dt            <- fChargingPumHydro -
        #dH_dt            <- #StoredHydro
        dS_dt            <- fSPCharging
        
        return(list(c(dL_dt, dP_dt, dS_dt),
                    
                    ChargingLithium = fChargingLithium,
                    ChargingPumpedHydro = fChargingPumHydro,
                    ExcessPower = fSPCharging
                    ))
      })
    }
    

    
    o <- data.frame(ode(y=stocks, times=simtime, func = model,
                        parms = auxs, method = "euler"))
    
    
    #Basically Plotting the data
    flow_plot <- ggplot(data = o, mapping = aes(time, ChargingLithium)) + theme_classic() +
      geom_line(data = o, mapping = aes(time, ChargingLithium), size = 1, color = "blue", linetype =2)+
      geom_line(data = o, mapping = aes(time, ChargingPumpedHydro), size = 1, color = "red",linetype =2)+
      geom_line(data = o, mapping = aes(time, ExcessPower), size = 1, color = "black")
    
    f <-   renderPlot({
      flow_plot <- ggplot(data = o, mapping = aes(time, ChargingLithium)) + theme_classic() +
        geom_line(data = o, mapping = aes(time, ChargingLithium), size = 1, color = "blue", linetype =2)+
        geom_line(data = o, mapping = aes(time, ChargingPumpedHydro), size = 1, color = "red",linetype =2)+
        geom_line(data = o, mapping = aes(time, ExcessPower), size = 1, color = "black")
    })
    
    capital_plot <- ggplot(data = o, mapping = aes(time, sLithium)) + theme_classic() +
      geom_line(data = o, mapping = aes(time, sLithium), size = 1, color = "blue", linetype =2)+
      geom_line(data = o, mapping = aes(time, sPumpedHydro), size = 1, color = "black")+
      geom_line(data = o, mapping = aes(time, sSuperPower), size = 1, color = "yellow")
    
    ressource_plot <- ggplot(data = o, mapping = aes(time, sCapital)) + theme_classic() +
      geom_line(data = o, mapping = aes(time, sResource), size = 1, color = "black", linetype =1)
    
    output$arrange <- renderPlot({
      grid.arrange(flow_plot,capital_plot,ressource_plot, nrow = 3)
    })
    
  })
}

shinyApp(ui, server)

Now the error I get is: Warning: Error in currentLithiumChargePower: object 'aCLithium' not found.

I am really not sure if I put the helper functions in the right spot or if I have to put them somewhere outside the observe function. Altough I think this is not the problem, as I tried with a print statement and the previous function chargeFunction does execute. The problem seems to be in the assignment:

    currentLithiumChargePower <- function(sLithium, aCharge, aLithiumCap, aEfficiencyLithium) {
  print("I execute!!!")
  availableCap <- (aLithiumCap-sLithium)
  potCharge    <- (aCharge/aEfficiencyLithium)
  chargePower  <- (aCLithium*aLithiumCap)/aEfficiencyLithium
  print("I don't execute")

I have the input data in a csv, but since it doesn't really matter I just generated random vectors. If you'd like the data though I'm happy to provide it.

Thank you very much for reading and hopefully you can help me out.

Can you please simplify your code, i.e. provide a **minimal** reproducible example, see https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example — tpetzoldt, Nov 04 '21 at 22:21

score 1 · Answer 1 · answered Nov 04 '21 at 22:27

The code provided by the OP is quite long, so we don't knof if someone invests time to it. In case the "helper function" is just a function that exists outside of server and ui, it can just be implemented as "global". Here a minimal example taken from here:

library("deSolve")

brusselator <- function(t, y, p) {
  with(as.list(c(y, p)), {
    dX <- k1*A   - k2*B*X    + k3*X^2*Y - k4*X
    dY <- k2*B*X - k3*X^2*Y
    list(c(X=dX, Y=dY))
  })
}

server <- function(input, output) {
  output$brussels <- renderPlot({
    parms <- c(A=input$A, B=input$B, k1=1, k2=1, k3=1, k4=1)
    out <- ode(y = c(X=1, Y=1), times=seq(0, 100, .1), brusselator, parms)
    matplot.0D(out)
  })
}

ui <- fluidPage(
  numericInput("A", label = "A", value = 1),
  numericInput("B", label = "B", value = 3),
  plotOutput("brussels")
)

shinyApp(ui=ui, server=server)

Note also that an observer function can be avoided here.

score 0 · Answer 2 · answered Nov 06 '21 at 21:32

Thanks a lot tpetzoldt for the answer.

I could solve my issue by just passing through the variable in the function call.

aLithiumCharging <- currentLithiumChargePower(sLithium,aCharge, aLithiumCap,aEfficiencyLithium,aCLithium)

Which I figured out after reducing my problem to the minimal replicable example. I thought I tried that but apparently not.

Thanks a lot

R [Shiny] How to implement helper functions into a dynamical model using deSolve and Shiny

2 Answers2