27 Jan 20:07
Unavailability modeling issue
John A Schroeder <John.Schroeder <at> inl.gov>
2010-01-27 19:07:28 GMT
2010-01-27 19:07:28 GMT
Hello All,
I am evaluating SimPy to see if it might be useful in my work. I have stumbled across something that looks like a software bug, but my just be my lack of understanding.
I am using a variation of one of the MachRep.py examples to test my understanding of SimPy. The issue I have been trying to resolve is how to deal with the situation where the last scheduled machine "upTime" or "downTime" is longer than the machine mission time (or simulation end time). When this occurs, as it will in almost every simulation, it will cause an error in the unavailability calculation. I have tried to limit the scheduled time using code like the following:
if (now() + self.downTime) > GlobalVars.simEndTime:
self.downTime = GlobalVars.simEndTime - now() - .0001
The thing I do not understand is why, without subtracting the .0001 in the above code, SimPy appears to enter some sort of infinite loop. I have appended the full source code from my test problem. If I remove the .0001 in the above example, the calculation appears to repeat the above statement forever. What is happening?
John A. Schroeder
Idaho National Laboratory
Battelle Energy Alliance, LLC
P.O. Box 1625
Idaho Falls, ID 83415-3850
Ph: 208-526-8755
FAX: 208-526-2930
# EPowerSim01.py ===============================================================
#
# One machine, which sometimes break down.
# Up time is exponentially distributed
# Down time is exponentially distributed
#
# Output is long-run fraction of down time (unavailability).
from SimPy.Simulation import Monitor, Process, hold, initialize
from SimPy.Simulation import activate, simulate, now
from SimPy.SimPlot import SimPlot
from random import Random
from math import exp
## Global data -----------------------------------------------------------------
class GlobalVars: # Global variables
failureRate = 1.0E-3 # machine failure rate
dtMon = Monitor('Down Times')
utMon = Monitor('Up Times')
uaMon = Monitor('Unavailability')
mUA = 0 # machine unavailability
repairRate = 1.0E-1 # machine repair rate
rnd = Random(12345) # random number stream
simEndTime = 8760.0 # simulation end time
trials = 5000 # number of simulations to perform
## Model components ------------------------------------------------------------
# Calculate the true unavailability of a machine with failure rate l, repair
# rate m, and mission time t.
def UA(l,m,t):
return (l /(l + m)) * (1.0 - exp(-1.0*(l + m) * t))
class Machine(Process):
def __init__(self, id, fr, rr): # Required constructor
Process.__init__(self) # must call parent constructor
# Instance variables
self.ID = id # ID for this Machine object
self.failureRate = fr # failure rate
self.repairRate = rr # repair rate
self.totalDownTime = 0 # Keep track of total down time
self.totalUpTime = 0 # Keep track of total up time
self.upTime = 0
self.downTime = 0
def getDT(self): # Return the down time
return self.totalDownTime
def getUA(self): # Return the machine unavailability
return self.totalDownTime / (self.totalUpTime + self.totalDownTime)
def getUT(self): # Return the total up time
return self.totalUpTime
def run(self): # Required Process Execution Method (PEM)
while True:
# Machine is in the "up" state, hold for
# exponentially distributed up time.
self.upTime = GlobalVars.rnd.expovariate(self.failureRate)
if (now() + self.upTime) > GlobalVars.simEndTime:
self.upTime = GlobalVars.simEndTime - now()
GlobalVars.utMon.observe(self.upTime)
# keep track of total up time
self.totalUpTime += self.upTime
yield hold,self,self.upTime
# Machine is now in the "down" state, hold for
# exponentially distributed repair time.
#print "%7.3E Machine %d failed" % (now(), self.ID)
self.downTime = GlobalVars.rnd.expovariate(self.repairRate)
#print " Down time pick: %7.3E" % (downTime)
if (now() + self.downTime) > GlobalVars.simEndTime:
self.downTime = GlobalVars.simEndTime - now() - .0001
#print " Now + down time: %7.3E" % (now() + self.downTime)
GlobalVars.dtMon.observe(self.downTime)
# keep track of total down time
self.totalDownTime += self.downTime
yield hold,self,self.downTime
#print "%7.3E Machine %d repaired " % (now(),self.ID)
## Experiment data -------------------------------------------------------------
## Model/Experiment ------------------------------------------------------------
def model(trial):
print "Starting trial", trial+1
initialize()
# Create the Machine object
edg1 = Machine(1, GlobalVars.failureRate, GlobalVars.repairRate)
# Register the machine thread and execute its run() method,
activate(edg1,edg1.run())
# Run until simulation reaches end time.
simulate(until=GlobalVars.simEndTime)
ut = edg1.getUT()
dt = edg1.getDT()
ua = dt / (ut+dt)
GlobalVars.uaMon.observe(ua)
print " Total up time: %7.3E" % (ut)
print " Total down time: %7.3E" % (dt)
print " Sum of up + down: %7.3E" % (ut+dt)
print " Unavailability: %7.3E" % (ua)
# Save the machine unavailability for late use.
GlobalVars.mUA = edg1.getUA()
print " "
def main():
for trial in range(GlobalVars.trials):
model(trial)
## Results ---------------------------------------------------------------------
print "Machine %d mean up time: %7.3E" % (1, GlobalVars.utMon.mean())
print "Machine %d mean down time: %7.3E" % (1, GlobalVars.dtMon.mean())
print "Machine %d mean UA: %7.3E" % (1, GlobalVars.uaMon.mean())
print "Machine %d true UA: %7.3E" % (1, UA(GlobalVars.failureRate,
GlobalVars.repairRate,
GlobalVars.simEndTime))
#histo = GlobalVars.dtMon.histogram(low=0.0, high=100. ,nbins=50)
#plt = SimPlot()
#plt.plotHistogram(histo,xlab='Time (hr)', ylab='Count',
# title="Down times",
# color="red",width=2)
#plt.mainloop()
if __name__ == '__main__': main()
------------------------------------------------------------------------------ The Planet: dedicated and managed hosting, cloud storage, colocation Stay online with enterprise data centers and the best network in the business Choose flexible plans and management services without long-term contracts Personal 24x7 support from experience hosting pros just a phone call away. http://p.sf.net/sfu/theplanet-com
_______________________________________________ Simpy-users mailing list Simpy-users <at> lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/simpy-users
RSS Feed