library(poweRlaw)
library(MonteCarlo)
#t<-read.csv("auto sim setup.txt")
setup<-read.csv("wbs_fit setup.csv")
nr=length(setup[[1]])
t=setup
ttf <- function(i) {
# i is machine index
v=NULL
#for (i in 1:nr) {
if ((t[i,6]=="D L-N")|(t[i,6]=="C L-N")|(t[i,6]=="fit L-N")) {
x<-rlnorm(1000,t[i,4],t[i,5])+t[i,3]
}
if ((t[i,6]=="C Wei")|(t[i,6]=="fit Wei")) {
x<-rweibull(1000,t[i,4],t[i,5])+t[i,3]
}
if ((t[i,6]=="D Exp")|(t[i,6]=="C Exp")|(t[i,6]=="fit Exp")) {
x<-rexp(1000,t[i,4])+t[i,3]
}
if (t[i,6]=="D Poi") {
x<-rpois(1000,t[i,4])+t[i,3]
}
if (is.null(v)) v=data.frame(x)
names(v)=c("v1")
return (list("mean"=mean(v$v1)))
}
param_list=list("i"=c(1:nr))
erg<-MonteCarlo(func=ttf, nrep=2500, param_list=param_list, ncpus=1)
## Grid of 6 parameter constellations to be evaluated.
##
## Progress:
##
##
|
| | 0%
## Warning in searchCommandline(parallel, cpus = cpus, type = type,
## socketHosts = socketHosts, : Unknown option on commandline:
## rmarkdown::render('C:/Users/itais/Documents/1R/dT with WBS/dt wbs predict
## faults for 2017&2018 v2.R', encoding
##
|
|=========== | 17%
## Warning in searchCommandline(parallel, cpus = cpus, type = type,
## socketHosts = socketHosts, : Unknown option on commandline:
## rmarkdown::render('C:/Users/itais/Documents/1R/dT with WBS/dt wbs predict
## faults for 2017&2018 v2.R', encoding
##
|
|====================== | 33%
## Warning in searchCommandline(parallel, cpus = cpus, type = type,
## socketHosts = socketHosts, : Unknown option on commandline:
## rmarkdown::render('C:/Users/itais/Documents/1R/dT with WBS/dt wbs predict
## faults for 2017&2018 v2.R', encoding
##
|
|================================ | 50%
## Warning in searchCommandline(parallel, cpus = cpus, type = type,
## socketHosts = socketHosts, : Unknown option on commandline:
## rmarkdown::render('C:/Users/itais/Documents/1R/dT with WBS/dt wbs predict
## faults for 2017&2018 v2.R', encoding
##
|
|=========================================== | 67%
## Warning in searchCommandline(parallel, cpus = cpus, type = type,
## socketHosts = socketHosts, : Unknown option on commandline:
## rmarkdown::render('C:/Users/itais/Documents/1R/dT with WBS/dt wbs predict
## faults for 2017&2018 v2.R', encoding
##
|
|====================================================== | 83%
## Warning in searchCommandline(parallel, cpus = cpus, type = type,
## socketHosts = socketHosts, : Unknown option on commandline:
## rmarkdown::render('C:/Users/itais/Documents/1R/dT with WBS/dt wbs predict
## faults for 2017&2018 v2.R', encoding
##
|
|=================================================================| 100%
##
df<-MakeFrame(erg)
renew_rate.mean=NULL
renew_rate.sd=NULL
f=file(description = "dt wbs ttf x_results.txt", open = "wt", encoding = "UTF-8")
for (i in 1:nr) {
x = df$mean[df$i==i]
renew_rate.mean[i]<-mean(x)
renew_rate.sd[i]<-sd(x)
hist(x = x,freq = F, xlab="TTF [days]", main=t[i,1])
legend(x = "bottom", legend = c("empirical (from data)","Normal distribution"), pch = c(0,0), col = c("black","red") )
lines(density(x), lty="dashed")
xnorm = rnorm(n = length(x), mean = renew_rate.mean[i], sd = renew_rate.sd[i])
lines(density(xnorm), col="red")
ks_result = numeric(2500)
for (j in 1:2500) {
w3 = ks.test(xnorm,x) #using KS test to see how good the fit is
if (w3$p.value > 0.10){
ks_result[j] = 1}
if (w3$p.value <= 0.10){
ks_result[j] = 0}
xnorm = rnorm(n = length(x), mean = renew_rate.mean[i], sd = renew_rate.sd[i])
}
print(paste0("Confidence of ",setup[i,1]," KS test:",sum(ks_result)/25,"%"))
if (sum(ks_result)>0.9*2500) { print("failed to reject the null hypothesis. renew rate is normally dist.") }
else { print("rejected the null hypothesis. renew rate is not normally dist.")}
write(paste(setup[i,1],",",renew_rate.mean[i],",",renew_rate.sd[i],"KS confidence,",sum(ks_result)/25,"%"),file = f)
}
## [1] "Confidence of Eagle.Cont KS test:96.48%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Eagle.Mech KS test:96.68%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Eagle.Elec KS test:97.64%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matrix.Cont KS test:97.44%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matrix.Mech KS test:97.04%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matrix.Elec KS test:98.04%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
close(f)
count.N<-function(mu, t) {
# mu is the mean time between Events
# t is the time period
# returns m(t) - mean number of events on [0,t]
n=0
s=mu # start position is the 1st event time
while (s<=t) {
n=n+1
s=s+mu
}
return (n)
}
# cl is a list of c which is the cost ratio of PM / ER
# or which is the ttr ratio of PM / ER
# rw is a list of PM periods in days
# eff is the machine efficiency
cl=seq(0.1,0.9,0.1)
#rw=c(20,30,40,50,60,70,80,90,120)
rw=c(36, 38, 59, 83, 142, 397)
ER=rep(0,nr)
eff=0.74
AV<-rep(0,nr)
AO<-rep(0,nr)
#### Q of eta.A
mn=matrix(rep(0,9*nr),9,nr,dimnames = list(cl,setup[,1])) # minimum cost of c per PM period
Q=matrix(rep(0,9*nr),9,nr,dimnames = list(cl,setup[,1])) # Q factor
tmn=matrix(rep(0,9*nr),9,nr,dimnames = list(cl,setup[,1])) # the PM period which derive minimum cost
AVmx=matrix(rep(0,9*nr),9,nr,dimnames = list(cl,setup[,1])) # maximum availability of c per PM period
tAVmx=matrix(rep(0,9*nr),9,nr,dimnames = list(cl,setup[,1])) # the PM period which derive maximum availability
# eta.A is the average cost per Cost-Type-Criterion
# eta.B is the average cost per Availability-Criterion
for (i in 1:nr) { # Machine/Responsibility
eta.A<-matrix(rep(0,9*length(rw)),length(rw),9,dimnames = list(rw,cl))
N.A<-matrix(rep(0,9*length(rw)),length(rw),9,dimnames = list(rw,cl))
ER[i]<-1/renew_rate.mean[i]
eta.B<-matrix(rep(0,9*length(rw)),length(rw),9,dimnames = list(rw,cl))
AV[i]<-(renew_rate.mean[i])/(renew_rate.mean[i]+(setup[i,7]/(24*eff)))
AO[i]<-(renew_rate.mean[i])/(renew_rate.mean[i]+(setup[i,7]/(24*eff))+setup[i,8])
for (c in cl) {
for (t in rw) {
eta.A[which(rw==t),which(cl==c)]=(c+count.N(renew_rate.mean[i],t))/t
N.A[which(rw==t),which(cl==c)]=count.N(renew_rate.mean[i],t)
eta.B[which(rw==t),which(cl==c)] = 1 / (1 + eta.A[which(rw==t),which(cl==c)]*(setup[i,7]/(24*eff)) )
}
}
plot(cl,eta.A[1,], xlim=c(0,1),ylim=c(0,max(eta.A)), type='b',pch=1,xlog=T, main=paste0(setup[i,1]," Mean cost per unit time"))
legend(x="topleft",legend = c(rw),col=c(1:length(rw)),pch=c(1:length(rw)))
for(ix in 2:length(rw)) {lines(cl,eta.A[ix,],col=ix,type='b', xlog=T, pch=ix)}
lines(cl,rep(ER[i],length(cl)),col=1,type='l',xlog=T,lty='dashed')
for (c in 1:9) { # PM/ER cost ratio
mn[c,i]=min(eta.A[,c])
for (j in 1:length(rw)) { # PM period
if (mn[c,i]==eta.A[j,c]) {
tmn[c,i]=rw[j]
break
}
}
Q[c,i]=ER[i]/mn[c,i]
}
plot(cl,eta.B[1,], xlim=c(0,1),ylim=c(min(eta.B,AO[i]),max(eta.B)), type='b',pch=1,xlog=T, main=paste0(setup[i,1]," Mean availability per unit time"))
legend(x="topleft",legend = c(rw),col=c(1:length(rw)),pch=c(1:length(rw)))
for(ix in 2:length(rw)) {lines(cl,eta.B[ix,],col=ix,type='b', xlog=T, pch=ix)}
lines(cl,rep(AV[i],length(cl)),col='black',type='l',xlog=T,lty='dashed')
lines(cl,rep(AO[i],length(cl)),col='red',type='l',xlog=T,lty='dashed')
for (c in 1:9) { # PM/ER cost ratio
AVmx[c,i]=max(eta.B[,c])
for (j in 1:length(rw)) { # PM period
if (AVmx[c,i]==eta.B[j,c]) {
tAVmx[c,i]=rw[j]
break
}
}
}
}
#### PM perfiod optimization technique:
#rw=sort(as.factor(0.65/ER[1:nr]))
#rw=as.numeric(as.character(rw))
#rw=as.numeric(format(rw,digits = 0))
### now add margin to rw by add 0.5*rw[1] in front of rw, and 1.5*rw[last] at the end
### run the loop again and check the impact of Q and tmn and mn
# see "finding best PM period.R"