library(fitdistrplus)
## Loading required package: MASS
## Loading required package: survival
library(FAdist)
library(MonteCarlo)
source("cost_functions.R")
## Loading required package: stats4
## Loading required package: splines
##
## Attaching package: 'VGAM'
## The following objects are masked from 'package:FAdist':
##
## dgev, dgumbel, pgev, pgumbel, qgev, qgumbel, rgev, rgumbel
setup<-read.csv("wbs_setup fitdistrplus weibull.csv")
setup=setup[!is.na(setup$type),]
nr=length(setup[[1]])
t=setup
ttf <- function(i) {
# i is machine index
v=NULL
mn=0
vr=0
if (t[i,3]!=0) {
x<-rweibull3(1000,t[i,4],t[i,5],t[i,3])
fit=fitdist(data = x,distr = "weibull")$estimate
mn=mean.wei(fit[1],fit[2])
vr=var.wei(fit[1],fit[2])
} else {
x<-rweibull(1000,t[i,4],t[i,5])
fit=fitdist(data = x,distr = "weibull")$estimate
mn=mean.wei(fit[1],fit[2])
vr=var.wei(fit[1],fit[2])
}
if (is.null(v)) v=data.frame(mn,sqrt(vr) )
names(v)=c("mean","sd")
return (list(mean=v$mean,sd=v$sd))
}
param_list=list("i"=c(1:nr))
erg<-MonteCarlo(func=ttf, nrep=500, param_list=param_list, ncpus=1)
## Grid of 34 parameter constellations to be evaluated.
##
## Progress:
##
##
|
| | 0%
|
|== | 3%
|
|==== | 6%
|
|====== | 9%
|
|======== | 12%
|
|========== | 15%
|
|=========== | 18%
|
|============= | 21%
|
|=============== | 24%
|
|================= | 26%
|
|=================== | 29%
|
|===================== | 32%
|
|======================= | 35%
|
|========================= | 38%
|
|=========================== | 41%
|
|============================= | 44%
|
|=============================== | 47%
|
|================================ | 50%
|
|================================== | 53%
|
|==================================== | 56%
|
|====================================== | 59%
|
|======================================== | 62%
|
|========================================== | 65%
|
|============================================ | 68%
|
|============================================== | 71%
|
|================================================ | 74%
|
|================================================== | 76%
|
|==================================================== | 79%
|
|====================================================== | 82%
|
|======================================================= | 85%
|
|========================================================= | 88%
|
|=========================================================== | 91%
|
|============================================================= | 94%
|
|=============================================================== | 97%
|
|=================================================================| 100%
##
df<-MakeFrame(erg)
renew_rate.mean=NULL
renew_rate.sd=NULL
f=file(description = "dt wbs ttf v4.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) #mean(df$sd[df$i==i])
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:97.96%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Eagle.Mech KS test:96.28%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matrix1.Cont KS test:98.12%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matrix1.Elec KS test:96.16%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matrix1.Gnrl KS test:97.64%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matrix1.Mech KS test:92.52%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matrix1.Smith KS test:95.8%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matrix2.Elec KS test:98.32%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matrix2.Mech KS test:98.6%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matrix2.Smith KS test:98.08%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto1.Cont KS test:85.24%"
## [1] "rejected the null hypothesis. renew rate is not normally dist."
## [1] "Confidence of Auto1.Elec KS test:98.16%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto1.Mech KS test:97.28%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto2.Cont KS test:96.92%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto2.Elec KS test:92.36%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto2.Mech KS test:98.28%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto2.Smith KS test:98.56%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto3.Cont KS test:95.72%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto3.Elec KS test:96.4%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto3.Mech KS test:98.28%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto4.Cont KS test:97.92%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto4.Elec KS test:96.88%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Auto4.Mech KS test:98.28%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Bullmer.Cont KS test:97.88%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Bullmer.Elec KS test:94.64%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Bullmer.Mech KS test:90.56%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of C-scan1.Elec KS test:98.56%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of C-scan1.Mech KS test:95.12%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of C-scan1.Smith KS test:95.44%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Eastman.Cont KS test:97.6%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Eastman.Mech KS test:97.6%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matec.Cont KS test:98.52%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matec.Elec KS test:98.88%"
## [1] "failed to reject the null hypothesis. renew rate is normally dist."
## [1] "Confidence of Matec.Cont KS test:97.84%"
## [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.2,1.0,0.2)
cln=length(cl)
#rw=c(20,30,40,50,60,70,80,90,120)
rw.mean=c(8,16,20,24,29,39,43,57,66)
rw.mean_and_sd=sort(c(372,35,134,56,79))
rw = rw.mean
ER=rep(0,nr)
ER.L=ER
ER.H=ER
eff=0.74
AV<-rep(0,nr)
AV.L=AV
AO<-rep(0,nr)
#### Q of eta.A
mn=matrix(rep(0,cln*nr),cln,nr,dimnames = list(cl,setup[,1])) # minimum cost of c per PM period
Q=matrix(rep(0,cln*nr),cln,nr,dimnames = list(cl,setup[,1])) # Q factor
tmn=matrix(rep(0,cln*nr),cln,nr,dimnames = list(cl,setup[,1])) # the PM period which derive minimum cost
AVmx=matrix(rep(0,cln*nr),cln,nr,dimnames = list(cl,setup[,1])) # maximum availability of c per PM period
tAVmx=matrix(rep(0,cln*nr),cln,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
span=2 #1 or 2 or 3 times the std dev
for (i in 1:nr) { # Machine/Responsibility
if (is.na(setup[i,7])) next
if (setup[i,2]!="IFR" & setup[i,2]!="EXP") next
eta.A<-matrix(rep(0,cln*length(rw)),length(rw),cln,dimnames = list(rw,cl))
N.A<-matrix(rep(0,cln*length(rw)),length(rw),cln,dimnames = list(rw,cl))
ER.L[i]<-1/(renew_rate.mean[i]-span*renew_rate.sd[i])
ER.H[i]<-1/(renew_rate.mean[i]+span*renew_rate.sd[i])
ER[i]<-1/(renew_rate.mean[i])
eta.B<-matrix(rep(0,cln*length(rw)),length(rw),cln,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])
AV.L[i]<-(renew_rate.mean[i]+span*renew_rate.sd[i])/(renew_rate.mean[i]+span*renew_rate.sd[i]+((setup[i,7]+setup[i,9])/(24*eff)))
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)) )
}
}
#png(filename = paste0(setup[i,1]," Mean cost per unit time.png"), width = 800, height = 600)
plot(cl,eta.A[1,], xlim=c(0,1),ylim=c(0,max(eta.A)), type='b',pch=1,
main=paste0(setup[i,1]," ",setup[i,2]," Mean cost per unit time"),
xlab="PM/ER cost ratio", ylab="Cost per unit time [d]")
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,ER[i]*cl,col='black',type='l',xlog=T,lty='dashed')
lines(cl,ER.L[i]*rep(1,cln),col='grey',type='l',xlog=T,lty='dashed')
lines(cl,ER.H[i]*rep(1,cln),col='grey',type='l',xlog=T,lty='dashed')
#dev.off()
for (c in 1:cln) { # 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]
}
#png(filename = paste0(setup[i,1]," Mean availability per unit time.png"), width = 800, height = 600)
plot(cl,eta.B[1,], xlim=c(0,1), ylim=c(min(eta.B,AO[i],AV.L[i]),max(eta.B)), type='b',pch=1,
main=paste0(setup[i,1]," ",setup[i,2]," Mean availability per unit time"),
xlab="PM/ER cost ratio",ylab="Availability")
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(AV.L[i],length(cl)),col='grey',type='l',xlog=T,lty='dashed')
lines(cl,rep(AO[i],length(cl)),col='red',type='l',xlog=T,lty='dashed')
#dev.off()
for (c in 1:cln) { # 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"