-- 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.Specs import Simulation.Aivika.Simulation import Simulation.Aivika.Dynamics import Simulation.Aivika.Event import Simulation.Aivika.Ref import Simulation.Aivika.Process import Simulation.Aivika.Random import qualified Simulation.Aivika.DoubleLinkedList as DLL import Simulation.Aivika.Experiment import Simulation.Aivika.Experiment.ExperimentSpecsView import Simulation.Aivika.Experiment.FinalStatsView import Simulation.Aivika.Experiment.DeviationChartView import Simulation.Aivika.Experiment.FinalHistogramView -- | 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 { furnaceNormalGen :: IO Double, -- ^ The normal random number generator. furnacePits :: [Pit], -- ^ The pits for ingots. furnacePitCount :: Ref Int, -- ^ The count of active pits with ingots. furnaceAwaitingIngots :: DLL.DoubleLinkedList Ingot, -- ^ The awaiting ingots in the queue. furnaceQueueCount :: Ref Int, -- ^ 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 { pitIngot :: Ref (Maybe Ingot), -- ^ The ingot in the pit. pitTemp :: Ref 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 Furnace newFurnace = do normalGen <- liftIO newNormalGen pits <- sequence [newPit | i <- [1..10]] pitCount <- newRef 0 awaitingIngots <- liftIO DLL.newList queueCount <- newRef 0 waitCount <- newRef 0 waitTime <- newRef 0.0 heatingTime <- newRef 0.0 h <- newRef 1650.0 totalCount <- newRef 0 loadCount <- newRef 0 unloadCount <- newRef 0 unloadTemps <- newRef [] return Furnace { furnaceNormalGen = normalGen, furnacePits = pits, furnacePitCount = pitCount, furnaceAwaitingIngots = awaitingIngots, furnaceQueueCount = queueCount, furnaceWaitCount = waitCount, furnaceWaitTime = waitTime, furnaceHeatingTime = heatingTime, furnaceTemp = h, furnaceTotalCount = totalCount, furnaceLoadCount = loadCount, furnaceUnloadCount = unloadCount, furnaceUnloadTemps = unloadTemps } -- | Create a new pit. newPit :: Simulation Pit newPit = do ingot <- newRef Nothing h' <- newRef 0.0 return Pit { pitIngot = ingot, pitTemp = h' } -- | Create a new ingot. newIngot :: Furnace -> Event Ingot newIngot furnace = do t <- liftDynamics 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 -> Event () 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' <- liftDynamics 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 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 () 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 -> Event () tryLoadPit furnace pit = do let ingots = furnaceAwaitingIngots furnace flag <- liftIO $ DLL.listNull ingots unless flag $ do ingot <- liftIO $ DLL.listFirst ingots liftIO $ DLL.listRemoveFirst ingots t' <- liftDynamics time modifyRef (furnaceQueueCount furnace) (+ (-1)) loadIngot (ingot { ingotLoadTime = t', ingotLoadTemp = 400.0 }) pit -- | Unload the ingot from the specified pit. unloadIngot :: Ingot -> Pit -> Event () 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 <- readRef (furnacePitCount furnace) writeRef (furnacePitCount furnace) (count - 1) -- how long did we heat the ingot up? t' <- liftDynamics 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 -> Event () loadIngot ingot pit = do writeRef (pitIngot pit) $ Just ingot writeRef (pitTemp pit) $ ingotLoadTemp ingot -- count the active pits let furnace = ingotFurnace ingot count <- readRef (furnacePitCount furnace) writeRef (furnacePitCount furnace) (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' <- liftDynamics time 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 -> Event () startIteratingFurnace furnace = let pits = furnacePits furnace in enqueueEventWithIntegTimes $ 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' <- liftDynamics dt h <- readRef (furnaceTemp furnace) writeRef (furnaceTemp furnace) $ h + dt' * (2600.0 - h) * 0.2 -- | Return all empty pits. emptyPits :: Furnace -> Event [Pit] emptyPits furnace = filterM (fmap isNothing . readRef . pitIngot) $ furnacePits furnace -- | Accept a new ingot. acceptIngot :: Furnace -> Event () acceptIngot furnace = do ingot <- newIngot furnace -- counting modifyRef (furnaceTotalCount furnace) (+ 1) -- check what to do with the new ingot count <- readRef (furnacePitCount furnace) if count >= 10 then do let ingots = furnaceAwaitingIngots furnace liftIO $ DLL.listAddLast ingots ingot modifyRef (furnaceQueueCount furnace) (+ 1) 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 liftEvent $ acceptIngot furnace -- repeat it again processFurnace furnace -- | Initialize the furnace. initializeFurnace :: Furnace -> Event () 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 furnace <- newFurnace pid <- newProcessId -- initialize the furnace and start its iterating in start time runEventInStartTime IncludingCurrentEvents $ do initializeFurnace furnace startIteratingFurnace furnace -- accept input ingots runProcessInStartTime IncludingCurrentEvents pid (processFurnace furnace) -- the mean wait time in the queue let meanWaitTime :: Event Double meanWaitTime = do waitTime <- readRef (furnaceWaitTime furnace) waitCount <- readRef (furnaceWaitCount furnace) return $ waitTime / fromIntegral waitCount -- the mean heating time let meanHeatingTime :: Event Double meanHeatingTime = do heatingTime <- readRef (furnaceHeatingTime furnace) unloadCount <- readRef (furnaceUnloadCount furnace) return $ heatingTime / fromIntegral unloadCount experimentDataInStartTime [(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), (pitCountName, seriesEntity "the used pit count" $ furnacePitCount furnace), (queueCountName, seriesEntity "the queue size" $ furnaceQueueCount furnace), (meanWaitTimeName, seriesEntity "the mean wait time" meanWaitTime), (meanHeatingTimeName, seriesEntity "the mean heating time" meanHeatingTime) ] totalIngotCountName = "totalIngotCount" loadedIngotCountName = "loadedIngotCount" readyIngotCountName = "readyIngotCount" awaitedIngotCountName = "awaitedIngotCount" readyIngotTempsName = "readyIngotTemps" pitCountName = "pitCount" queueCountName = "queueCount" 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 pitCountName] }, outputView $ defaultFinalHistogramView { finalHistogramTitle = "Final Histogram - 2", finalHistogramPlotTitle = "The distribution of the used pit count " ++ "in the final simulation time point.", finalHistogramSeries = [pitCountName] }, outputView $ defaultFinalStatsView { finalStatsTitle = "Final Statistics - 2", finalStatsDescription = "The summary of the used pit count in the final time point.", finalStatsSeries = [pitCountName] }, outputView $ defaultDeviationChartView { deviationChartTitle = "Deviation Chart - 3", deviationChartPlotTitle = "The queue size", deviationChartSeries = [Right queueCountName] }, outputView $ defaultFinalHistogramView { finalHistogramTitle = "Final Histogram - 3", finalHistogramPlotTitle = "The distribution of the queue size " ++ "in the final simulation time point.", finalHistogramSeries = [queueCountName] }, outputView $ defaultFinalStatsView { finalStatsTitle = "Final Statistics - 3", finalStatsDescription = "The summary of the queue size in the final time point.", finalStatsSeries = [queueCountName] }, 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