diff --git a/Simulation/Aivika/Branch/Generator.hs b/Simulation/Aivika/Branch/Generator.hs --- a/Simulation/Aivika/Branch/Generator.hs +++ b/Simulation/Aivika/Branch/Generator.hs @@ -22,6 +22,7 @@ import Data.IORef import Simulation.Aivika.Trans +import Simulation.Aivika.Trans.Generator.Primitive import Simulation.Aivika.Branch.Internal.Br instance MonadGenerator BrIO where @@ -37,8 +38,12 @@ generateUniformInt = generateUniformInt01 . generator01 + generateTriangular = generateTriangular01 . generator01 + generateNormal = generateNormal01 . generatorNormal01 + generateLogNormal = generateLogNormal01 . generatorNormal01 + generateExponential = generateExponential01 . generator01 generateErlang = generateErlang01 . generator01 @@ -47,6 +52,14 @@ generateBinomial = generateBinomial01 . generator01 + generateGamma g = generateGamma01 (generatorNormal01 g) (generator01 g) + + generateBeta g = generateBeta01 (generatorNormal01 g) (generator01 g) + + generateWeibull = generateWeibull01 . generator01 + + generateDiscrete = generateDiscrete01 . generator01 + newGenerator tp = case tp of SimpleGenerator -> @@ -71,44 +84,6 @@ return Generator { generator01 = g01, generatorNormal01 = gNormal01 } --- | Generate an uniform random number with the specified minimum and maximum. -generateUniform01 :: BrIO Double - -- ^ the generator - -> Double - -- ^ minimum - -> Double - -- ^ maximum - -> BrIO Double -generateUniform01 g min max = - do x <- g - return $ min + x * (max - min) - --- | Generate an uniform random number with the specified minimum and maximum. -generateUniformInt01 :: BrIO Double - -- ^ the generator - -> Int - -- ^ minimum - -> Int - -- ^ maximum - -> BrIO Int -generateUniformInt01 g min max = - do x <- g - let min' = fromIntegral min - max' = fromIntegral max - return $ round (min' + x * (max' - min')) - --- | Generate a normal random number by the specified generator, mean and variance. -generateNormal01 :: BrIO Double - -- ^ normal random numbers with mean 0 and variance 1 - -> Double - -- ^ mean - -> Double - -- ^ variance - -> BrIO Double -generateNormal01 g mu nu = - do x <- g - return $ mu + nu * x - -- | Create a normal random number generator with mean 0 and variance 1 -- by the specified generator of uniform random numbers from 0 to 1. newNormalGenerator01 :: BrIO Double @@ -148,63 +123,3 @@ liftIO $ writeIORef flagRef True liftIO $ writeIORef nextRef $ xi2 * psi return $ xi1 * psi - --- | Return the exponential random number with the specified mean. -generateExponential01 :: BrIO Double - -- ^ the generator - -> Double - -- ^ the mean - -> BrIO Double -generateExponential01 g mu = - do x <- g - return (- log x * mu) - --- | Return the Erlang random number. -generateErlang01 :: BrIO Double - -- ^ the generator - -> Double - -- ^ the scale - -> Int - -- ^ the shape - -> BrIO Double -generateErlang01 g beta m = - do x <- loop m 1 - return (- log x * beta) - where loop m acc - | m < 0 = error "Negative shape: generateErlang." - | m == 0 = return acc - | otherwise = do x <- g - loop (m - 1) (x * acc) - --- | Generate the Poisson random number with the specified mean. -generatePoisson01 :: BrIO Double - -- ^ the generator - -> Double - -- ^ the mean - -> BrIO Int -generatePoisson01 g mu = - do prob0 <- g - let loop prob prod acc - | prob <= prod = return acc - | otherwise = loop - (prob - prod) - (prod * mu / fromIntegral (acc + 1)) - (acc + 1) - loop prob0 (exp (- mu)) 0 - --- | Generate a binomial random number with the specified probability and number of trials. -generateBinomial01 :: BrIO Double - -- ^ the generator - -> Double - -- ^ the probability - -> Int - -- ^ the number of trials - -> BrIO Int -generateBinomial01 g prob trials = loop trials 0 where - loop n acc - | n < 0 = error "Negative number of trials: generateBinomial." - | n == 0 = return acc - | otherwise = do x <- g - if x <= prob - then loop (n - 1) (acc + 1) - else loop (n - 1) acc diff --git a/Simulation/Aivika/Branch/Internal/Br.hs b/Simulation/Aivika/Branch/Internal/Br.hs --- a/Simulation/Aivika/Branch/Internal/Br.hs +++ b/Simulation/Aivika/Branch/Internal/Br.hs @@ -1,4 +1,6 @@ +{-# LANGUAGE RecursiveDo #-} + -- | -- Module : Simulation.Aivika.Branch.Internal.Br -- Copyright : Copyright (c) 2016, David Sorokin @@ -24,6 +26,7 @@ import Control.Applicative import Control.Monad import Control.Monad.Trans +import Control.Monad.Fix import Control.Exception (throw, catch, finally) import Simulation.Aivika.Trans.Exception @@ -76,6 +79,12 @@ {-# INLINE liftIO #-} liftIO = Br . const . liftIO + +instance MonadFix BrIO where + + mfix f = + Br $ \ps -> + do { rec { a <- invokeBr ps (f a) }; return a } instance MonadException BrIO where diff --git a/aivika-branches.cabal b/aivika-branches.cabal --- a/aivika-branches.cabal +++ b/aivika-branches.cabal @@ -1,5 +1,5 @@ name: aivika-branches -version: 0.1 +version: 0.1.1 synopsis: Branching discrete event simulation library description: This package extends the Aivika [1] library with facilities for creating branches to run @@ -19,6 +19,8 @@ tested-with: GHC == 7.10.3 extra-source-files: tests/MachRep1.hs + tests/Furnace.hs + tests/InspectionAdjustmentStations.hs library @@ -37,7 +39,7 @@ containers >= 0.4.0.0, random >= 1.0.0.3, aivika >= 4.3.2, - aivika-transformers >= 4.3.1 + aivika-transformers >= 4.3.2 extensions: TypeFamilies, MultiParamTypeClasses, diff --git a/tests/Furnace.hs b/tests/Furnace.hs new file mode 100644 --- /dev/null +++ b/tests/Furnace.hs @@ -0,0 +1,325 @@ + +-- This is a model of the Furnace. It is described in different sources [1, 2]. +-- +-- [1] A. Alan B. Pritsker, Simulation with Visual SLAM and AweSim, 2nd ed. +-- +-- [2] Труб И.И., Объектно-ориентированное моделирование на C++: Учебный курс. - СПб.: Питер, 2006 + +import Data.Maybe +import Control.Monad +import Control.Monad.Trans + +import Simulation.Aivika.Trans +import Simulation.Aivika.Trans.Queue.Infinite +import Simulation.Aivika.Branch + +type DES = BrIO + +-- | The simulation specs. +specs = Specs { spcStartTime = 0.0, + -- spcStopTime = 1000.0, + spcStopTime = 300.0, + spcDT = 0.1, + spcMethod = RungeKutta4, + spcGeneratorType = SimpleGenerator } + +-- | Return a random initial temperature of the item. +randomTemp :: Parameter DES Double +randomTemp = randomUniform 400 600 + +-- | Represents the furnace. +data Furnace = + Furnace { furnacePits :: [Pit], + -- ^ The pits for ingots. + furnacePitCount :: Ref DES Int, + -- ^ The count of active pits with ingots. + furnaceQueue :: FCFSQueue DES Ingot, + -- ^ The furnace queue. + furnaceUnloadedSource :: SignalSource DES (), + -- ^ Notifies when the ingots have been + -- unloaded from the furnace. + furnaceHeatingTime :: Ref DES (SamplingStats Double), + -- ^ The heating time for the ready ingots. + furnaceTemp :: Ref DES Double, + -- ^ The furnace temperature. + furnaceReadyCount :: Ref DES Int, + -- ^ The count of ready ingots. + furnaceReadyTemps :: Ref DES [Double] + -- ^ The temperatures of all ready ingots. + } + +-- | Notifies when the ingots have been unloaded from the furnace. +furnaceUnloaded :: Furnace -> Signal DES () +furnaceUnloaded = publishSignal . furnaceUnloadedSource + +-- | A pit in the furnace to place the ingots. +data Pit = + Pit { pitIngot :: Ref DES (Maybe Ingot), + -- ^ The ingot in the pit. + pitTemp :: Ref DES Double + -- ^ The ingot temperature in the pit. + } + +data Ingot = + Ingot { ingotFurnace :: Furnace, + -- ^ The furnace. + ingotReceiveTime :: Double, + -- ^ The time at which the ingot was received. + ingotReceiveTemp :: Double, + -- ^ The temperature with which the ingot was received. + ingotLoadTime :: Double, + -- ^ The time of loading in the furnace. + ingotLoadTemp :: Double, + -- ^ The temperature when the ingot was loaded in the furnace. + ingotCoeff :: Double + -- ^ The heating coefficient. + } + +-- | Create a furnace. +newFurnace :: Simulation DES Furnace +newFurnace = + do pits <- sequence [newPit | i <- [1..10]] + pitCount <- newRef 0 + queue <- runEventInStartTime newFCFSQueue + heatingTime <- newRef emptySamplingStats + h <- newRef 1650.0 + readyCount <- newRef 0 + readyTemps <- newRef [] + s <- newSignalSource + return Furnace { furnacePits = pits, + furnacePitCount = pitCount, + furnaceQueue = queue, + furnaceUnloadedSource = s, + furnaceHeatingTime = heatingTime, + furnaceTemp = h, + furnaceReadyCount = readyCount, + furnaceReadyTemps = readyTemps } + +-- | Create a new pit. +newPit :: Simulation DES Pit +newPit = + do ingot <- newRef Nothing + h' <- newRef 0.0 + return Pit { pitIngot = ingot, + pitTemp = h' } + +-- | Create a new ingot. +newIngot :: Furnace -> Event DES Ingot +newIngot furnace = + do t <- liftDynamics time + xi <- liftParameter $ randomNormal 0.05 0.01 + h' <- liftParameter randomTemp + let c = 0.1 + xi + return Ingot { ingotFurnace = furnace, + ingotReceiveTime = t, + ingotReceiveTemp = h', + ingotLoadTime = t, + ingotLoadTemp = h', + ingotCoeff = c } + +-- | Heat the ingot up in the pit if there is such an ingot. +heatPitUp :: Pit -> Event DES () +heatPitUp pit = + do ingot <- readRef (pitIngot pit) + case ingot of + Nothing -> + return () + Just ingot -> do + + -- update the temperature of the ingot. + let furnace = ingotFurnace ingot + dt' <- liftParameter dt + h' <- readRef (pitTemp pit) + h <- readRef (furnaceTemp furnace) + writeRef (pitTemp pit) $ + h' + dt' * (h - h') * ingotCoeff ingot + +-- | Check whether there are ready ingots in the pits. +ingotsReady :: Furnace -> Event DES Bool +ingotsReady furnace = + fmap (not . null) $ + filterM (fmap (>= 2200.0) . readRef . pitTemp) $ + furnacePits furnace + +-- | Try to unload the ready ingot from the specified pit. +tryUnloadPit :: Furnace -> Pit -> Event DES () +tryUnloadPit furnace pit = + do h' <- readRef (pitTemp pit) + when (h' >= 2000.0) $ + do Just ingot <- readRef (pitIngot pit) + unloadIngot furnace ingot pit + +-- | Try to load an awaiting ingot in the specified empty pit. +tryLoadPit :: Furnace -> Pit -> Event DES () +tryLoadPit furnace pit = + do ingot <- tryDequeue (furnaceQueue furnace) + case ingot of + Nothing -> + return () + Just ingot -> + do t' <- liftDynamics time + loadIngot furnace (ingot { ingotLoadTime = t', + ingotLoadTemp = 400.0 }) pit + +-- | Unload the ingot from the specified pit. +unloadIngot :: Furnace -> Ingot -> Pit -> Event DES () +unloadIngot furnace ingot pit = + do h' <- readRef (pitTemp pit) + writeRef (pitIngot pit) Nothing + writeRef (pitTemp pit) 0.0 + + -- count the active pits + modifyRef (furnacePitCount furnace) (+ (- 1)) + + -- how long did we heat the ingot up? + t' <- liftDynamics time + modifyRef (furnaceHeatingTime furnace) $ + addSamplingStats (t' - ingotLoadTime ingot) + + -- what is the temperature of the unloaded ingot? + modifyRef (furnaceReadyTemps furnace) (h' :) + + -- count the ready ingots + modifyRef (furnaceReadyCount furnace) (+ 1) + +-- | Load the ingot in the specified pit +loadIngot :: Furnace -> Ingot -> Pit -> Event DES () +loadIngot furnace ingot pit = + do writeRef (pitIngot pit) $ Just ingot + writeRef (pitTemp pit) $ ingotLoadTemp ingot + + -- count the active pits + modifyRef (furnacePitCount furnace) (+ 1) + count <- readRef (furnacePitCount furnace) + + -- decrease the furnace temperature + h <- readRef (furnaceTemp furnace) + let h' = ingotLoadTemp ingot + dh = - (h - h') / fromIntegral count + writeRef (furnaceTemp furnace) $ h + dh + +-- | Start iterating the furnace processing through the event queue. +startIteratingFurnace :: Furnace -> Event DES () +startIteratingFurnace furnace = + let pits = furnacePits furnace + in enqueueEventWithIntegTimes $ + do -- try to unload ready ingots + ready <- ingotsReady furnace + when ready $ + do mapM_ (tryUnloadPit furnace) pits + triggerSignal (furnaceUnloadedSource furnace) () + + -- heat up + mapM_ heatPitUp pits + + -- update the temperature of the furnace + dt' <- liftParameter dt + h <- readRef (furnaceTemp furnace) + writeRef (furnaceTemp furnace) $ + h + dt' * (2600.0 - h) * 0.2 + +-- | Return all empty pits. +emptyPits :: Furnace -> Event DES [Pit] +emptyPits furnace = + filterM (fmap isNothing . readRef . pitIngot) $ + furnacePits furnace + +-- | This process takes ingots from the queue and then +-- loads them in the furnace. +loadingProcess :: Furnace -> Process DES () +loadingProcess furnace = + do ingot <- dequeue (furnaceQueue furnace) + let wait = + do count <- liftEvent $ readRef (furnacePitCount furnace) + when (count >= 10) $ + do processAwait (furnaceUnloaded furnace) + wait + wait + -- take any empty pit and load it + liftEvent $ + do pit: _ <- emptyPits furnace + loadIngot furnace ingot pit + -- repeat it again + loadingProcess furnace + +-- | The input process that adds new ingots to the queue. +inputProcess :: Furnace -> Process DES () +inputProcess furnace = + do delay <- liftParameter $ + randomExponential 2.5 + holdProcess delay + -- we have got a new ingot + liftEvent $ + do ingot <- newIngot furnace + enqueue (furnaceQueue furnace) ingot + -- repeat it again + inputProcess furnace + +-- | Initialize the furnace. +initializeFurnace :: Furnace -> Event DES () +initializeFurnace furnace = + do x1 <- newIngot furnace + x2 <- newIngot furnace + x3 <- newIngot furnace + x4 <- newIngot furnace + x5 <- newIngot furnace + x6 <- newIngot furnace + let p1 : p2 : p3 : p4 : p5 : p6 : ps = + furnacePits furnace + loadIngot furnace (x1 { ingotLoadTemp = 550.0 }) p1 + loadIngot furnace (x2 { ingotLoadTemp = 600.0 }) p2 + loadIngot furnace (x3 { ingotLoadTemp = 650.0 }) p3 + loadIngot furnace (x4 { ingotLoadTemp = 700.0 }) p4 + loadIngot furnace (x5 { ingotLoadTemp = 750.0 }) p5 + loadIngot furnace (x6 { ingotLoadTemp = 800.0 }) p6 + writeRef (furnaceTemp furnace) 1650.0 + +-- | The simulation model. +model :: Simulation DES (Results DES) +model = + do furnace <- newFurnace + + -- initialize the furnace and start its iterating in start time + runEventInStartTime $ + do initializeFurnace furnace + startIteratingFurnace furnace + + -- generate randomly new input ingots + runProcessInStartTime $ + inputProcess furnace + + -- load permanently the input ingots in the furnace + runProcessInStartTime $ + loadingProcess furnace + + -- return the simulation results + return $ + resultSummary $ + results + [resultSource "inputIngotCount" "the input ingot count" $ + enqueueStoreCount (furnaceQueue furnace), + -- + resultSource "loadedIngotCount" "the loaded ingot count" $ + dequeueCount (furnaceQueue furnace), + -- + resultSource "outputIngotCount" "the output ingot count" $ + furnaceReadyCount furnace, + -- + resultSource "outputIngotTemp" "the output ingot temperature" $ + fmap listSamplingStats $ readRef $ furnaceReadyTemps furnace, + -- + resultSource "heatingTime" "the heating time" $ + furnaceHeatingTime furnace, + -- + resultSource "pitCount" "the number of ingots in pits" $ + furnacePitCount furnace, + -- + resultSource "furnaceQueue" "the furnace queue" $ + furnaceQueue furnace] + +-- | The main program. +main = + runBr $ + printSimulationResultsInStopTime + printResultSourceInEnglish + model specs diff --git a/tests/InspectionAdjustmentStations.hs b/tests/InspectionAdjustmentStations.hs new file mode 100644 --- /dev/null +++ b/tests/InspectionAdjustmentStations.hs @@ -0,0 +1,168 @@ + +{-# LANGUAGE RecursiveDo, Arrows #-} + +-- Example: Inspection and Adjustment Stations on a Production Line +-- +-- This is a model of the workflow with a loop. Also there are two infinite queues. +-- +-- It is described in different sources [1, 2]. So, this is chapter 8 of [2] and section 5.15 of [1]. +-- +-- [1] A. Alan B. Pritsker, Simulation with Visual SLAM and AweSim, 2nd ed. +-- +-- [2] Труб И.И., Объектно-ориентированное моделирование на C++: Учебный курс. - СПб.: Питер, 2006 + +import Prelude hiding (id, (.)) + +import Control.Monad +import Control.Monad.Trans +import Control.Monad.Fix +import Control.Arrow +import Control.Category (id, (.)) + +import Simulation.Aivika.Trans +import Simulation.Aivika.Trans.Queue.Infinite +import Simulation.Aivika.Branch + +type DES = BrIO + +-- | The simulation specs. +specs = Specs { spcStartTime = 0.0, + spcStopTime = 480.0, + spcDT = 0.1, + spcMethod = RungeKutta4, + spcGeneratorType = SimpleGenerator } + +-- the minimum delay of arriving the next TV set +minArrivalDelay = 3.5 + +-- the maximum delay of arriving the next TV set +maxArrivalDelay = 7.5 + +-- the minimum time to inspect the TV set +minInspectionTime = 6 + +-- the maximum time to inspect the TV set +maxInspectionTime = 12 + +-- the probability of passing the inspection phase +inspectionPassingProb = 0.85 + +-- how many are inspection stations? +inspectionStationCount = 2 + +-- the minimum time to adjust an improper TV set +minAdjustmentTime = 20 + +-- the maximum time to adjust an improper TV set +maxAdjustmentTime = 40 + +-- how many are adjustment stations? +adjustmentStationCount = 1 + +-- create an inspection station (server) +newInspectionStation :: Simulation DES (Server DES () a (Either a a)) +newInspectionStation = + newServer $ \a -> + do holdProcess =<< + (liftParameter $ + randomUniform minInspectionTime maxInspectionTime) + passed <- + liftParameter $ + randomTrue inspectionPassingProb + if passed + then return $ Right a + else return $ Left a + +-- create an adjustment station (server) +newAdjustmentStation :: Simulation DES (Server DES () a a) +newAdjustmentStation = + newServer $ \a -> + do holdProcess =<< + (liftParameter $ + randomUniform minAdjustmentTime maxAdjustmentTime) + return a + +model :: Simulation DES (Results DES) +model = mdo + -- to count the arrived TV sets for inspecting and adjusting + inputArrivalTimer <- newArrivalTimer + -- it will gather the statistics of the processing time + outputArrivalTimer <- newArrivalTimer + -- define a stream of input events + let inputStream = + randomUniformStream minArrivalDelay maxArrivalDelay + -- create a queue before the inspection stations + inspectionQueue <- + runEventInStartTime newFCFSQueue + -- create a queue before the adjustment stations + adjustmentQueue <- + runEventInStartTime newFCFSQueue + -- create the inspection stations (servers) + inspectionStations <- + forM [1 .. inspectionStationCount] $ \_ -> + newInspectionStation + -- create the adjustment stations (servers) + adjustmentStations <- + forM [1 .. adjustmentStationCount] $ \_ -> + newAdjustmentStation + -- a processor loop for the inspection stations' queue + let inspectionQueueProcessorLoop = + queueProcessorLoopSeq + (liftEvent . enqueue inspectionQueue) + (dequeue inspectionQueue) + inspectionProcessor + (adjustmentQueueProcessor >>> adjustmentProcessor) + -- a processor for the adjustment stations' queue + let adjustmentQueueProcessor = + queueProcessor + (liftEvent . enqueue adjustmentQueue) + (dequeue adjustmentQueue) + -- a parallel work of the inspection stations + let inspectionProcessor = + processorParallel (map serverProcessor inspectionStations) + -- a parallel work of the adjustment stations + let adjustmentProcessor = + processorParallel (map serverProcessor adjustmentStations) + -- the entire processor from input to output + let entireProcessor = + arrivalTimerProcessor inputArrivalTimer >>> + inspectionQueueProcessorLoop >>> + arrivalTimerProcessor outputArrivalTimer + -- start simulating the model + runProcessInStartTime $ + sinkStream $ runProcessor entireProcessor inputStream + -- return the simulation results in start time + return $ + results + [resultSource + "inspectionQueue" "the inspection queue" + inspectionQueue, + -- + resultSource + "adjustmentQueue" "the adjustment queue" + adjustmentQueue, + -- + resultSource + "inputArrivalTimer" "the input arrival timer" + inputArrivalTimer, + -- + resultSource + "outputArrivalTimer" "the output arrival timer" + outputArrivalTimer, + -- + resultSource + "inspectionStations" "the inspection stations" + inspectionStations, + -- + resultSource + "adjustmentStations" "the adjustment stations" + adjustmentStations] + +modelSummary :: Simulation DES (Results DES) +modelSummary = fmap resultSummary model + +main = + runBr $ + printSimulationResultsInStopTime + printResultSourceInEnglish + modelSummary specs