John A Schroeder | 27 Jan 20:07
Favicon

Unavailability modeling issue


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

Gmane