#Compute average year
datasim=scan("datasim.txt")
datasim=ts(datasim,frequency=365)
pr1=cycle(datasim)
avgyear=tapply(datasim,pr1,mean)

#Plot average year
par(mar=c(5,5,2,2),cex=1.5)
plot(avgyear,type="l",xlab="Days",ylab=expression('Average daily inflow (m'^3*'/s)'))
grid(col="black")

#Plot cumulative inflow volume
par(mar=c(5,5,2,2),cex=1.5)
plot(cumsum(avgyear)*86400/1000000,type="l",xlab="Days",ylab=expression('Cumulative inflow volume (10'^6*'m'^3*')'),lwd=2)
grid(col="black")

#Compute number of failures
vol=15000000
volt=rep(0,18250)
volt[1]=7000000
release=50
fail=0
for (i in 2:18250)
{
volt[i]=volt[i-1]+datasim[i]*86400-release*86400
if(volt[i]>vol)volt[i]=vol
if(volt[i]<0)
{
volt[i]=0
fail=fail+1
}
}

#Compute number of failures through a function
contaf=function(release)
{
volt=rep(0,18250)
fail=0
volt[1]=7000000
for (i in 2:18250)
{
volt[i]=volt[i-1]+datasim[i]*86400-release*86400
if(volt[i]>vol)volt[i]=vol
if(volt[i]<0)
{
volt[i]=0
fail=fail+1
}
}
return(fail)
}

#Compute the benefit
benefit=function(x,a=350,b=0.01,c=10,d=0.1,corr=1)
{
benefit=a*(1-exp(-b*x))-contaf(x)*corr/50*c*x^d
return(benefit)
}

#Maximum benefit should be obtained with a release between 220 and 230

#Maximisation of utility

ben1=0
for (i in 1:8)
ben1[i]=benefit(160+(i-1)*10)

ben2=0
for (i in 1:8)
ben2[i]=benefit(160+(i-1)*10,corr=0.5)

ben3=0
for (i in 1:8)
ben3[i]=benefit(160+(i-1)*10,corr=2)

mintot=min(ben1,ben2,ben3)

uben1=1-exp(-0.02*(ben1-mintot))
uben2=1-exp(-0.02*(ben2-mintot))
uben3=1-exp(-0.02*(ben3-mintot))

ures=0
for (i in 1:7)
ures[i]=mean(c(uben1[i],uben2[i],uben3[i]))

ures