-- This is a model of the Furnace. It is described in different sources [1, 2]. -- -- [1] { add a foreign source in English } -- -- [2] Труб И.И., Объектно-ориентированное моделирование на C++: Учебный курс. - СПб.: Питер, 2006 -- -- This model is often used in the literature as an example of combined -- continuous-discrete simulation but this is not a point here. It illustrates -- how the time-driven and process-oriented simulation models can be combined -- based on the common event queue. It still uses the differential equation but -- it is modeled directly [3] with help of the Euler method within the time-driven -- part of the combined model. -- -- [3] The time bounds for such an equation are much smaller than that ones which are defined -- by the specs. Therefore there is no sense to use the 'integ' function as it would be -- very slow because of large allocating memory for each integral, although it is possible. -- -- However, you can still combine the differential (and difference) equations with the DES and -- agent-based models. The integral (as well as any 'Dynamics' computation) can be used directly -- in the DES sub-model. But to update something from the DES sub-model that could be used aready -- in the differential equations, you should save data with help of types 'Var' or 'UVar' as they -- keep all the history of their past values. Also the values of these two types are managed by -- the event queue that allows synchronizing them with the DES sub-model. -- -- To define the external parameters for the Monte-Carlo simulation, see the Financial model. -- -- To enable the parallel simulation, you should compile it -- with option -threaded and then pass in other options +RTS -N2 -RTS -- to the executable if you have a dual core processor without -- hyper-threading. Also you can increase the number -- of parallel threads via option -N if you have a more modern -- processor. import Data.Maybe import System.Random import Control.Monad import Control.Monad.Trans import Simulation.Aivika.Dynamics import Simulation.Aivika.Dynamics.Simulation import Simulation.Aivika.Dynamics.Base import Simulation.Aivika.Dynamics.EventQueue import Simulation.Aivika.Dynamics.Ref import Simulation.Aivika.Dynamics.UVar import Simulation.Aivika.Dynamics.Process import Simulation.Aivika.Dynamics.Random import Simulation.Aivika.Statistics import Simulation.Aivika.Experiment import Simulation.Aivika.Experiment.ExperimentSpecsView import Simulation.Aivika.Experiment.FinalStatsView import Simulation.Aivika.Experiment.DeviationChartView import Simulation.Aivika.Experiment.FinalHistogramView import qualified Simulation.Aivika.Queue as Q -- | The simulation specs. specs = Specs { spcStartTime = 0.0, -- spcStopTime = 1000.0, spcStopTime = 300.0, spcDT = 0.1, spcMethod = RungeKutta4 } -- | Return an exponentially distributed random value with mean -- 1 / @lambda@, where @lambda@ is a parameter of the function. exprnd :: Double -> IO Double exprnd lambda = do x <- getStdRandom random return (- log x / lambda) -- | Return a random initial temperature of the item. temprnd :: IO Double temprnd = do x <- getStdRandom random return (400.0 + (600.0 - 400.0) * x) -- | Represents the furnace. data Furnace = Furnace { furnaceQueue :: EventQueue, -- ^ The event queue. furnaceNormalGen :: IO Double, -- ^ The normal random number generator. furnacePits :: [Pit], -- ^ The pits for ingots. furnacePitCount :: UVar Int, -- ^ The count of active pits with ingots. furnacePitCountStats :: Ref (SamplingStats Int), -- ^ The statistics about the active pits. furnaceAwaitingIngots :: Q.Queue Ingot, -- ^ The awaiting ingots in the queue. furnaceQueueCount :: UVar Int, -- ^ The queue count. furnaceQueueCountStats :: Ref (SamplingStats Int), -- ^ The statistics about the queue count. furnaceWaitCount :: Ref Int, -- ^ The count of awaiting ingots. furnaceWaitTime :: Ref Double, -- ^ The wait time for all loaded ingots. furnaceHeatingTime :: Ref Double, -- ^ The heating time for all unloaded ingots. furnaceTemp :: Ref Double, -- ^ The furnace temperature. furnaceTotalCount :: Ref Int, -- ^ The total count of ingots. furnaceLoadCount :: Ref Int, -- ^ The count of loaded ingots. furnaceUnloadCount :: Ref Int, -- ^ The count of unloaded ingots. furnaceUnloadTemps :: Ref [Double] -- ^ The temperatures of all unloaded ingots. } -- | A pit in the furnace to place the ingots. data Pit = Pit { pitQueue :: EventQueue, -- ^ The bound dynamics queue. pitIngot :: Ref (Maybe Ingot), -- ^ The ingot in the pit. pitTemp :: Ref Double -- ^ The ingot temperature in the pit. } data Ingot = Ingot { ingotFurnace :: Furnace, -- ^ Return 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 :: EventQueue -> Simulation Furnace newFurnace queue = do normalGen <- liftIO normalGen pits <- sequence [newPit queue | i <- [1..10]] pitCount <- newUVar queue 0 pitCountStats <- newRef queue emptySamplingStats awaitingIngots <- liftIO Q.newQueue queueCount <- newUVar queue 0 queueCountStats <- newRef queue emptySamplingStats waitCount <- newRef queue 0 waitTime <- newRef queue 0.0 heatingTime <- newRef queue 0.0 h <- newRef queue 1650.0 totalCount <- newRef queue 0 loadCount <- newRef queue 0 unloadCount <- newRef queue 0 unloadTemps <- newRef queue [] return Furnace { furnaceQueue = queue, furnaceNormalGen = normalGen, furnacePits = pits, furnacePitCount = pitCount, furnacePitCountStats = pitCountStats, furnaceAwaitingIngots = awaitingIngots, furnaceQueueCount = queueCount, furnaceQueueCountStats = queueCountStats, furnaceWaitCount = waitCount, furnaceWaitTime = waitTime, furnaceHeatingTime = heatingTime, furnaceTemp = h, furnaceTotalCount = totalCount, furnaceLoadCount = loadCount, furnaceUnloadCount = unloadCount, furnaceUnloadTemps = unloadTemps } -- | Create a new pit. newPit :: EventQueue -> Simulation Pit newPit queue = do ingot <- newRef queue Nothing h' <- newRef queue 0.0 return Pit { pitQueue = queue, pitIngot = ingot, pitTemp = h' } -- | Create a new ingot. newIngot :: Furnace -> Dynamics Ingot newIngot furnace = do t <- time xi <- liftIO $ furnaceNormalGen furnace h' <- liftIO temprnd let c = 0.1 + (0.05 + xi * 0.01) 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 -> Dynamics () 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' <- 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 -> Dynamics 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 -> Dynamics () tryUnloadPit furnace pit = do h' <- readRef (pitTemp pit) when (h' >= 2000.0) $ do Just ingot <- readRef (pitIngot pit) unloadIngot ingot pit -- | Try to load an awaiting ingot in the specified empty pit. tryLoadPit :: Furnace -> Pit -> Dynamics () tryLoadPit furnace pit = do let ingots = furnaceAwaitingIngots furnace flag <- liftIO $ Q.queueNull ingots unless flag $ do ingot <- liftIO $ Q.queueFront ingots liftIO $ Q.dequeue ingots t' <- time modifyUVar (furnaceQueueCount furnace) (+ (-1)) c <- readUVar (furnaceQueueCount furnace) modifyRef (furnaceQueueCountStats furnace) $ addSamplingStats c loadIngot (ingot { ingotLoadTime = t', ingotLoadTemp = 400.0 }) pit -- | Unload the ingot from the specified pit. unloadIngot :: Ingot -> Pit -> Dynamics () unloadIngot ingot pit = do h' <- readRef (pitTemp pit) writeRef (pitIngot pit) Nothing writeRef (pitTemp pit) 0.0 -- count the active pits let furnace = ingotFurnace ingot count <- readUVar (furnacePitCount furnace) writeUVar (furnacePitCount furnace) (count - 1) modifyRef (furnacePitCountStats furnace) $ addSamplingStats (count - 1) -- how long did we heat the ingot up? t' <- time modifyRef (furnaceHeatingTime furnace) (+ (t' - ingotLoadTime ingot)) -- what is the temperature of the unloaded ingot? modifyRef (furnaceUnloadTemps furnace) (h' :) -- count the unloaded ingots modifyRef (furnaceUnloadCount furnace) (+ 1) -- | Load the ingot in the specified pit loadIngot :: Ingot -> Pit -> Dynamics () loadIngot ingot pit = do writeRef (pitIngot pit) $ Just ingot writeRef (pitTemp pit) $ ingotLoadTemp ingot -- count the active pits let furnace = ingotFurnace ingot count <- readUVar (furnacePitCount furnace) writeUVar (furnacePitCount furnace) (count + 1) modifyRef (furnacePitCountStats furnace) $ addSamplingStats (count + 1) -- decrease the furnace temperature h <- readRef (furnaceTemp furnace) let h' = ingotLoadTemp ingot dh = - (h - h') / fromInteger (toInteger (count + 1)) writeRef (furnaceTemp furnace) $ h + dh -- how long did we keep the ingot in the queue? t' <- time when (ingotReceiveTime ingot < t') $ do modifyRef (furnaceWaitCount furnace) (+ 1) modifyRef (furnaceWaitTime furnace) (+ (t' - ingotReceiveTime ingot)) -- count the loaded ingots modifyRef (furnaceLoadCount furnace) (+ 1) -- | Start iterating the furnace processing through the event queue. startIteratingFurnace :: Furnace -> Dynamics () startIteratingFurnace furnace = let queue = furnaceQueue furnace pits = furnacePits furnace in enqueueWithIntegTimes queue $ do ready <- ingotsReady furnace when ready $ do mapM_ (tryUnloadPit furnace) pits pits' <- emptyPits furnace mapM_ (tryLoadPit furnace) pits' mapM_ heatPitUp pits -- update the temperature of the furnace dt' <- dt h <- readRef (furnaceTemp furnace) writeRef (furnaceTemp furnace) $ h + dt' * (2600.0 - h) * 0.2 -- | Return all empty pits. emptyPits :: Furnace -> Dynamics [Pit] emptyPits furnace = filterM (fmap isNothing . readRef . pitIngot) $ furnacePits furnace -- | Accept a new ingot. acceptIngot :: Furnace -> Dynamics () acceptIngot furnace = do ingot <- newIngot furnace -- counting modifyRef (furnaceTotalCount furnace) (+ 1) -- check what to do with the new ingot count <- readUVar (furnacePitCount furnace) if count >= 10 then do let ingots = furnaceAwaitingIngots furnace liftIO $ Q.enqueue ingots ingot modifyUVar (furnaceQueueCount furnace) (+ 1) c <- readUVar (furnaceQueueCount furnace) modifyRef (furnaceQueueCountStats furnace) $ addSamplingStats c else do pit:_ <- emptyPits furnace loadIngot ingot pit -- | Process the furnace. processFurnace :: Furnace -> Process () processFurnace furnace = do delay <- liftIO $ exprnd (1.0 / 2.5) holdProcess delay -- we have got a new ingot liftDynamics $ acceptIngot furnace -- repeat it again processFurnace furnace -- | Initialize the furnace. initializeFurnace :: Furnace -> Dynamics () 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 (x1 { ingotLoadTemp = 550.0 }) p1 loadIngot (x2 { ingotLoadTemp = 600.0 }) p2 loadIngot (x3 { ingotLoadTemp = 650.0 }) p3 loadIngot (x4 { ingotLoadTemp = 700.0 }) p4 loadIngot (x5 { ingotLoadTemp = 750.0 }) p5 loadIngot (x6 { ingotLoadTemp = 800.0 }) p6 writeRef (furnaceTotalCount furnace) 6 writeRef (furnaceTemp furnace) 1650.0 -- | The simulation model that returns experimental data. model :: Simulation ExperimentData model = do queue <- newQueue furnace <- newFurnace queue pid <- newProcessID queue -- initialize the furnace and start its iterating in start time runDynamicsInStartTime $ do initializeFurnace furnace startIteratingFurnace furnace -- start accepting input ingots by launching the process runDynamicsInStartTime $ do t0 <- starttime runProcess (processFurnace furnace) pid t0 experimentDataInStartTime queue $ [(totalIngotCountName, seriesEntity "total ingot count" $ furnaceTotalCount furnace), (loadedIngotCountName, seriesEntity "loaded ingot count" $ furnaceLoadCount furnace), (readyIngotCountName, seriesEntity "ready ingot count" $ furnaceUnloadCount furnace), (awaitedIngotCountName, seriesEntity "awaited in the queue ingot count" $ furnaceWaitCount furnace), (readyIngotTempsName, seriesEntity "the temperature of ready ingot" $ furnaceUnloadTemps furnace), (pitCountStatsName, seriesEntity "the used pit count" $ furnacePitCountStats furnace), (queueCountStatsName, seriesEntity "the queue size" $ furnaceQueueCountStats furnace), (meanWaitTimeName, seriesEntity "the mean wait time" $ readRef (furnaceWaitTime furnace) / fmap fromIntegral (readRef (furnaceWaitCount furnace))), (meanHeatingTimeName, seriesEntity "the mean heating time" $ readRef (furnaceHeatingTime furnace) / fmap fromIntegral (readRef (furnaceUnloadCount furnace))) ] totalIngotCountName = "totalIngotCount" loadedIngotCountName = "loadedIngotCount" readyIngotCountName = "readyIngotCount" awaitedIngotCountName = "awaitedIngotCount" readyIngotTempsName = "readyIngotTemps" pitCountStatsName = "pitCountStats" queueCountStatsName = "queueCountStats" meanWaitTimeName = "the mean wait time in the queue" meanHeatingTimeName = "the mean heating time" -- | The experiment. experiment :: Experiment experiment = defaultExperiment { experimentSpecs = specs, -- experimentRunCount = 1000, experimentRunCount = 100, experimentTitle = "The Furnace model (the Monte-Carlo simulation)", experimentGenerators = [outputView defaultExperimentSpecsView, outputView $ defaultDeviationChartView { deviationChartTitle = "Deviation Chart - 1.1", deviationChartPlotTitle = "The total, loaded and ready ingot counts", deviationChartSeries = [Right totalIngotCountName, Right loadedIngotCountName, Right readyIngotCountName] }, outputView $ defaultDeviationChartView { deviationChartTitle = "Deviation Chart - 1.2", deviationChartPlotTitle = "The awaited in the queue ingot count", deviationChartSeries = [Right awaitedIngotCountName] }, outputView $ defaultFinalHistogramView { finalHistogramTitle = "Final Histogram - 1.1", finalHistogramPlotTitle = "The distribution of total, loaded and ready " ++ "ingot counts in the final time point.", finalHistogramSeries = [totalIngotCountName, loadedIngotCountName, readyIngotCountName] }, outputView $ defaultFinalHistogramView { finalHistogramTitle = "Final Histogram - 1.2", finalHistogramPlotTitle = "The distribution of the awaited in the queue " ++ "ingot count in the final time point.", finalHistogramSeries = [awaitedIngotCountName] }, outputView $ defaultFinalStatsView { finalStatsTitle = "Final Statistics - 1", finalStatsDescription = "The summary of total, loaded, ready and awaited in " ++ "the queue ingot counts in the final time point.", finalStatsSeries = [totalIngotCountName, loadedIngotCountName, readyIngotCountName, awaitedIngotCountName] }, outputView $ defaultDeviationChartView { deviationChartTitle = "Deviation Chart - 2", deviationChartPlotTitle = "The used pit count", deviationChartSeries = [Right pitCountStatsName] }, outputView $ defaultFinalStatsView { finalStatsTitle = "Final Statistics - 2", finalStatsDescription = "The summary of the used pit count in the final time point.", finalStatsSeries = [pitCountStatsName] }, outputView $ defaultDeviationChartView { deviationChartTitle = "Deviation Chart - 3", deviationChartPlotTitle = "The queue size", deviationChartSeries = [Right queueCountStatsName] }, outputView $ defaultFinalStatsView { finalStatsTitle = "Final Statistics - 3", finalStatsDescription = "The summary of the queue size in the final time point.", finalStatsSeries = [queueCountStatsName] }, outputView $ defaultDeviationChartView { deviationChartTitle = "Deviation Chart - 4", deviationChartPlotTitle = "The mean wait time", deviationChartSeries = [Right meanWaitTimeName] }, outputView $ defaultFinalHistogramView { finalHistogramTitle = "Final Histogram - 4", finalHistogramPlotTitle = "The distribution of the mean wait time " ++ "in the final simulation time point.", finalHistogramSeries = [meanWaitTimeName] }, outputView $ defaultFinalStatsView { finalStatsTitle = "Final Statistics - 4", finalStatsDescription = "The summary of the mean wait time in " ++ "the final simulation time point.", finalStatsSeries = [meanWaitTimeName] }, outputView $ defaultDeviationChartView { deviationChartTitle = "Deviation Chart - 5", deviationChartPlotTitle = "The mean heating time", deviationChartSeries = [Right meanHeatingTimeName] }, outputView $ defaultFinalHistogramView { finalHistogramTitle = "Final Histogram - 5", finalHistogramPlotTitle = "The distribution of the mean heating time " ++ "in the final simulation time point.", finalHistogramSeries = [meanHeatingTimeName] }, outputView $ defaultFinalStatsView { finalStatsTitle = "Final Statistics - 5", finalStatsDescription = "The summary of the mean heating time in " ++ "the final simulation time point.", finalStatsSeries = [meanHeatingTimeName] }, outputView $ defaultDeviationChartView { deviationChartTitle = "Deviation Chart - 6", deviationChartPlotTitle = "The ready ingot temperature", deviationChartSeries = [Right readyIngotTempsName] }, outputView $ defaultFinalHistogramView { finalHistogramTitle = "Final Histogram - 6", finalHistogramPlotTitle = "The distribution of the ready ingot temperature " ++ "in the final simulation time point.", finalHistogramSeries = [readyIngotTempsName] }, outputView $ defaultFinalStatsView { finalStatsTitle = "Final Statistics - 6", finalStatsDescription = "The summary of the ready ingot temperature in " ++ "the final simulation time point.", finalStatsSeries = [readyIngotTempsName] } ] } -- | The main program that launches the simulation experiment to produce the HTML file. main = runExperimentParallel experiment model