library(ggplot2) library(dplyr) library(deSolve) library(gridExtra) library(cowplot) library(lubridate) #signal <- data.frame(times = times, import = rep(0, length(times))) Text=read.csv("input.csv") #head(signal) #signal$import <- ifelse((trunc(signal$times) %% 2 == 0), 0, 1) #signal[8:12,] #head(signal) df <- read.csv("data2020.csv",col.names=c("DateTime", "DB", "RH","DP","P","WIND","WIND-DIR","SUN","RAIN","START")) head(df) df=df[complete.cases(df), ] df$DateTime=as.POSIXct(df$DateTime,format="%d/%m/%Y %H:%M") df$DB=(df$DB/10) df <- df %>% filter(hour(DateTime) %in% c(0, 6, 12, 18),minute(DateTime) == 0) head(df,n=20) #stop() input <- approxfun(df$DB, rule = 2) model <- function(t, y, parms) { with(as.list(c(parms, y)), { Text <- input(t) # <---- here SP=20 HRS=ifelse(Text>y[1], max(0,min(0.7,(Text-y[1])*0.5)), max(0,min(0.7,(20-y[1])*1))) # HRS=ifelse(Text<16,0.7,0) dy1 <- -((qv * pa * Ca * (y[1] - ((HRS*y[1])-(Text*(HRS-1))))) / (V * pa * Ca))+((ifelse(Text>20,80,0)*(30*60))/(V*Ca*pa)) #+min(0,min(1,(y[1]-24)*1))*((q2 * pa * Ca * (9-y[1])) / (V * pa * Ca)) # dy2 = 1 # dy2 <- y[2] list(c(dy1), Text = Text, ifelse(Text>y[1], max(0,min(1,(Text-y[1])*0.5)), max(0,(min(1,(20-y[1])*1)))), ifelse(Text>y[1],0.5,max(0,min(1,(y[1]-20)*1)))*((q2 * pa * Ca * (9-y[1]))),qv*pa*Ca*(20-Text)) }) } parms = c(qv = 0.005*(30*60), pa = 1.204, Ca = 1006, V = 3, q2=0.005) out <- ode(y = c(Tr=20), times = seq(0,nrow(df),1), func = model, parms,rtol = 1e-1) out=as.data.frame(out) head(out,n=50) #parms = c(qv = 54, pa = 1.204, Ca = 1006, V = 3000) #model <- function(time, Tr,Text, parms) { #df <- read.csv("dataWEEK.csv",col.names=c("DateTime", "DB", "RH","DP","P","WIND","WIND-DIR","SUN","RAIN","START")) #df=df[complete.cases(df), ] #df$DateTime=as.Date(df$DateTime) #df$DB=(df$DB/10) # #Text=parms[5] #dTr <- -(parms[1] * parms[2] * parms[3] * (Tr - Text)) / (parms[4] * parms[2] * parms[3]) #dText = df$DB #list(dTr,Text) #}# #solution <- ode(y = 20, times =as.numeric(rownames(df)), func = model, parms = parms) #print(solution) #solution=as.data.frame(solution) names(out) <- c("time", "Tr","Text","Sig","Heating","Loss") out$Total <- cumsum((out$Heating/2)) head(out) ##head(solution) # # Sample data for plot1 # Create ggplot objects plot1 <- ggplot(out, aes(x=time, y=Tr)) + geom_line(color="red")+ geom_line(data=out, aes(x=time,y=Text),color="blue")+ #geom_hline(yintercept = 18)+ geom_hline(yintercept = 20) plot2 = ggplot(out, aes(x=time,y=Sig))+geom_bar(stat = "identity", width = 0.1, color = "black") plot3 = ggplot(out, aes(x=time,y=Heating))+geom_line()+ geom_line(data=out,aes(x=time,y=Loss),color="red") plot4 <- ggplot(out, aes(x=time, y=Total)) + geom_line() #plot5 = ggplot(out, aes(x=time, y=Error))+geom_line() # Open a PDF device to save the plots pdf("tx5.pdf", width = 40, height = 20) # Arrange and print the plots using gridExtra grid.arrange(plot1,plot2,plot3,plot4,nrow = 4) # Or arrange and print the plots using cowplot # plot_grid(plot1, plot2, nrow = 1) # Close the PDF device to save the plots to the PDF file dev.off() #gg_plot <- ggplot(data = out, aes(x = time, y = Total)) + # geom_line() + # geom_line(aes(x = time, y = Text), color = "blue")+ # geom_line(aes(x=time,y=Sig),color="red")+ # labs(x = "Time Steps", y = "Room Temperature (Tr)") + # ggtitle("Room Temperature Over Time") ## #ggsave(filename = "room_temperature_plot.pdf", plot = gg_plot)