R [ Shiny ]如何使用deSolve和Shiny将辅助函数实现到动态模型中

Question

我正在尝试实现我使用R-Shiny和deSolve包构建的Vensim模型。该模型本身的灵感来自Rethink X的一份报告，但我试图将其调整为适合瑞士的情况。

它基本上模拟了一个完全由风能、太阳能、锂电池和水力发电供电的能源电网。

我在R-Shiny中运行模型时遇到了麻烦，因为我从来没有在R中使用过helper函数和动态模型。

代码如下：

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)

现在我得到的错误是: Warning: error in currentLithiumChargePower:找不到对象'aCLithium‘。

我真的不确定我是把助手函数放在正确的位置，还是必须把它们放在观察函数之外的某个地方。尽管我认为这不是问题所在，因为我尝试使用了一条print语句，并且执行了前面的函数chargeFunction。问题似乎出在作业中：

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

我有一个csv格式的输入数据，但因为这并不重要，所以我只是生成了随机向量。如果你想要这些数据，我很乐意提供。

非常感谢你的阅读，希望你能帮助我。

Answer 1

OP提供的代码相当长，所以我们不知道是否有人在这上面投入了时间。如果“帮助函数”只是一个存在于server和ui之外的函数，则可以将其实现为“全局”。下面是取自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)

还要注意，这里可以避免使用observer函数。

Answer 2

非常感谢tpetzoldt的回答。

我可以通过在函数调用中传递变量来解决我的问题。

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

这是我在将问题简化为最小的可复制示例后得出的结论。我以为我试过了，但显然没有。

非常感谢