#lang scribble/doc
@(require
scribble/manual
scribblings/icons
(for-label racket/base
racket/contract
plot
(planet williams/simulation/simulation-with-graphics)))
@title[#:tag "continuous-simulation"]{Continuous Simulation Models}
@local-table-of-contents[]
The simulation collection supports continuous simulation models as well as discrete event models. Unlike discrete event simulations in which time advances in response to the times of events, continuous simulations advance time continuously by integrating the differential equations defining the system of equations implementing the continuous model.
The simulation collection uses the ordinary differential equation (ODE) solver from the science collection to implement its integration loop for continuous models. Many aspects of the integration loop may be controlled by the model developer, including the ODE step-type, ODE control type, and step size limitations.
A continuous model is represented by a process which has a @scheme[wait/continuously] statement that specifies the differential equations for the continuous model and, optionally, the termination condition for the model.
@section{Continuous Variables}
Continuous variables are used to define the data for continuous simulation models. Each continuous variable is associated with the process instance that created it and when that process instance executes a @scheme[work/continuously] statement, its continuous variables are added to the state vector maintained by the integrator.
Continuous variables have the same automatic data collection available as normal variables. That is, you may specify automatic data collection of continuous variables.
@defproc[(make-continuous-variable (initial-value real? 'uninitialized))
variable?]{
Returns a new continuous variable with the specified @scheme[initial-value]. If no @scheme[initial-value] is specified, @scheme['uninitialized] is used. The continuous variable is associated with the current process instance. A @scheme[make-continuous-variable] call is invalid outside of a process.}
The variable field @schemefont{continuous?} specifies whether or not a variable is a continuous variable.
The variable field @schemefont{state-index} is only used for continuous variables. When the process instance owning the continuous variable is in the @schemefont{PROCESS-WORKING-CONTINUOUSLY} state (that is, is currently in a @scheme[work/continuously]), this field contains the index of the variable in the state vector. In this state, the value of the variable and its derivative are maintained by the integration loop and variable values (using @scheme[get-variable-value]) are retrieved from the state vector. The @schemefont{state-index} field is @math{-1} for discrete variables or when a continuous variable is not being controlled by the integration loop.
@defproc[(variable-dt! (variable variable?)) real?]{
Returns the value of the derivative of @scheme[variable]. This value is retrieved from the derivative vector. An error is raised if the variable is not being controlled by the integration loop (that is, if the process owning the variable is not in the @schemefont{PROCESS-WORKING-CONTINUOUSLY} state).}
@defproc[(set-variable-dt! (variable variable?) (value real?)) any]{
Sets the value of the derivative of @scheme[variable] to @scheme[value]. The value is stored in the derivative vector. An error is raised if the variable is not being controlled by the integration loop (that is, if the process owning the variable is not in the @schemefont{PROCESS-WORKING-CONTINUOUSLY} state). This should be used in the @scheme[work/continuously] to implement the differential equations for the continuous model.}
@section{Simulation Control (Continuous)}
The simulation collection main loop (i.e., the @scheme[start-simulation] function) implements the integration loop internally.
A continuous simulation model is implemented within a process using the @scheme[work/continuously] statement. This is analogous to the @scheme[wait/work] statement in a discrete simulation model, but specifies the differential equations and, optionally, the terminating condition for the continuous model.
@defform*[#:id work/continuously #:literals (until)
((work/continuously
(until expr)
body ...)
(work/continuously
body ...))]{
Defines a continuous model within a process instance. The @scheme[body] expressions should compute the values of the derivatives of the relevant continuous variables in the process and set them using the @scheme[set-variable-dt!] function. Note that side effects in the @scheme[body] expressions should be avoided since the integration loop will typically evaluate the @scheme[body] expressions multiple times for each time step and may even regress in time when the error estimates in the integration loop are too high.
The @scheme[expr] in the @scheme[until] clause is evaluated at the start of each time step. If the @scheme[expr] is true, the process drops out of the @scheme[work/continuously] and resumes execution as a discrete event simulation model.}
@section{Example---Furnace Model 2}
This combined discrete-continuous simulation model is derived from an example in @italic{Introduction to Combined Discrete-Continuous Simulation Using SIMSCRIPT II.5} by Abdel-Moaty M Fayek @cite["Fayek02"].
@schememod[
racket/base
(code:comment "Model 2 - Continuous Simulation Model")
(require (planet williams/simulation/simulation-with-graphics))
(code:comment "Simulation Parameters")
(define end-time 720.0)
(define n-pits 7)
(code:comment "Data collection variables")
(define total-ingots 0)
(define wait-time #f)
(define heat-time #f)
(define leave-temp #f)
(code:comment "Model Definition")
(define random-sources (make-random-source-vector 4))
(define furnace-set #f)
(define furnace-temp 1500.0)
(define pit #f)
(define (scheduler)
(for ((i (in-naturals)))
(schedule #:now (ingot i))
(wait (random-exponential (vector-ref random-sources 0) 1.5))))
(define-process (ingot i)
(let* ((initial-temp (random-flat (vector-ref random-sources 1) 100.0 200.0))
(heat-coeff (+ (random-gaussian
(vector-ref random-sources 2) 0.05 0.01) 0.07))
(final-temp (random-flat (vector-ref random-sources 3) 800.0 1000.0))
(current-temp (make-continuous-variable initial-temp))
(arrive-time (current-simulation-time))
(start-time #f))
(with-resource (pit)
(when (= (modulo i 100) 0)
(accumulate (variable-history current-temp)))
(set-variable-value! wait-time (- (current-simulation-time) arrive-time))
(set-insert! furnace-set self)
(set! start-time (current-simulation-time))
(work/continuously
until (>= (variable-value current-temp) final-temp)
(set-variable-dt!
current-temp
(* (- furnace-temp (variable-value current-temp))
heat-coeff)))
(set-variable-value! heat-time (- (current-simulation-time) start-time))
(set-variable-value! leave-temp (variable-value current-temp))
(set-remove! furnace-set self))
(when (variable-history current-temp)
(write-special (history-plot (variable-history current-temp)
(format "Ingot ~a Temp History" i)))
(newline))
(set! total-ingots (+ total-ingots 1))))
(define (stop-sim)
(printf "Report after ~a Simulated Hours - ~a Ingots Processed~n"
(current-simulation-time) total-ingots)
(printf "~n-- Ingot Waiting Time Statistics --~n")
(printf "Mean Wait Time = ~a~n" (variable-mean wait-time))
(printf "Variance = ~a~n" (variable-variance wait-time))
(printf "Maximum Wait Time = ~a~n" (variable-maximum wait-time))
(printf "~n-- Ingot Heating Time Statistics --~n")
(printf "Mean Heat Time = ~a~n" (variable-mean heat-time))
(printf "Variance = ~a~n" (variable-variance heat-time))
(printf "Maximum Heat Time = ~a~n" (variable-maximum heat-time))
(printf "Minimum Heat Time = ~a~n" (variable-minimum heat-time))
(printf "~n-- Final Temperature Statistics --~n")
(printf "Mean Leave Temp = ~a~n" (variable-mean leave-temp))
(printf "Variance = ~a~n" (variable-variance leave-temp))
(printf "Maximum Leave Temp = ~a~n" (variable-maximum leave-temp))
(printf "Minimum Leave Temp = ~a~n" (variable-minimum leave-temp))
(write-special (history-plot (variable-history leave-temp)
"Final Temperature Histogram"))
(newline)
(printf "~n-- Furnace Utilization Statistics --~n")
(printf "Mean No. of Ingots = ~a~n"
(variable-mean (set-variable-n furnace-set)))
(printf "Variance = ~a~n"
(variable-variance (set-variable-n furnace-set)))
(printf "Maximum No. of Ingots = ~a~n"
(variable-maximum (set-variable-n furnace-set)))
(printf "Minimum No. of Ingots = ~a~n"
(variable-minimum (set-variable-n furnace-set)))
(write-special (history-plot (variable-history (set-variable-n furnace-set))
"Furnace Utilization History"))
(newline)
(stop-simulation))
(define (initialize)
(set! total-ingots 0)
(set! wait-time (make-variable))
(set! heat-time (make-variable))
(set! leave-temp (make-variable))
(set! pit (make-resource n-pits))
(set! furnace-set (make-set))
(accumulate (variable-history (set-variable-n furnace-set)))
(tally (variable-statistics wait-time))
(tally (variable-statistics heat-time))
(tally (variable-statistics leave-temp))
(tally (variable-history leave-temp))
(schedule #:at 0.0 (scheduler))
(schedule #:at end-time (stop-sim)))
(define (run-simulation)
(with-new-simulation-environment
(initialize)
(start-simulation)))
(run-simulation)
]
The following is the output from the model.
@image["scribblings/images/model-2-plot-1.gif"]
...
@image["scribblings/images/model-2-plot-2.gif"]
@verbatim{
Report after 720.0 Simulated Hours - 479 Ingots Processed
-- Ingot Waiting Time Statistics --
Mean Wait Time = 0.44596828009621853
Variance = 1.004363939227031
Maximum Wait Time = 6.380898029191428
-- Ingot Heating Time Statistics --
Mean Heat Time = 7.032338489600859
Variance = 1.0511669721776045
Maximum Heat Time = 10.88442744599513
Minimum Heat Time = 4.867835442878146
-- Final Temperature Statistics --
Mean Leave Temp = 912.2921783644567
Variance = 3300.2864186657825
Maximum Leave Temp = 1020.1172721368804
Minimum Leave Temp = 804.2299580285976}
@image["scribblings/images/model-2-plot-3.gif"]
@verbatim{
-- Furnace Utilization Statistics --
Mean No. of Ingots = 4.687638930905194
Variance = 3.321384522948101
Maximum No. of Ingots = 7
Minimum No. of Ingots = 0}
@image["scribblings/images/model-2-plot-4.gif"]
@section{Simulation Environment (Continuous)}
There are several field in each simulation environment that are specific to continuous models.
The @schemefont{continuous-event-list} field maintains a list of the current continuous events. Continuous event execution is controlled by the integration loop. Basically, the events on the continuous event list point to the body of the corresponding @scheme[work/continuously] statement. These (should) implement the continuous models by computing the appropriate variable derivatives. This field is set internally.
The @schemefont{evolve} field contains the ordinary differential equation evolver object. The evolver object combines the results of the stepping function and the control function, if specified, to reliably advance the simulation forward over the appropriate time interval. If the control function signals that the step size should be decreased, the evolution function backs out of the current step and tries the proposed smaller step size. This process is continued until an acceptable step size is found. This field is set internally.
The @schemefont{control} field contains the ordinary differential equation control object of false, @scheme[#f], for a fixed step size. The control function, if specified, examines the proposed changes to the solution and its error estimate produced by a stepping function and attempts to determine the optimal step size for a user-specified level of error. The default control function keeps the local error within an absolute error of @math{10}@superscript{@math{-6}} and relative error of @math{0.0} with respect to the solution. This field may be set by the user---in particular, set this field to false, @scheme[#f], for a fixed step size.
The @schemefont{step-type} field contains the ordinary differential equation solver step-type object. The step-type object specifies the stepping algorithm. The default step-type algorithm is an embedded (classical) Runge-Kutta-Fehberg (4, 5) method, which is a good general-purpose integrator. This field may be set by the user.
The @schemefont{step} field contains the ordinary equation stepper function. The stepper function advances the solution from time @math{t} to time @math{t + h} for a fixed step size @math{h} and estimates the resulting local error. This field is set internally.
The @schemefont{system} field contains the ordinary differential equation system of equations. This is created internally from the continuous variables in the processes currently working continuously (i.e., on the continuous event list). This field is set internally.
The @schemefont{step-size} field contains the current (last) step size for the ordinary differential equation solver. If the @schemefont{control} field is false, @scheme{#f}, this specifies the fixed step size; otherwise, it is controlled by the control object. The default step size is @math{10}@superscript{@math{-6}} (only applicable when the @schemefont{control} field is false, @scheme[#f]). This field may be set by the user.
The @schemefont{dimension} field contains the dimensionality of the system of equations for the ordinary differential solver. This is set internally.
The @schemefont{y} field contains the state vector for the system of equations for the ordinary differential solver. This is created internally from the continuous variables in the processes currently working continuously (i.e., on the continuous event list). The @scheme[get-variable-value] function retrieves the values of variables from the state vector when they are controlled by the integration loop. [Note that this field is false, @scheme[#f], when the integration loop is not running.]
The @schemefont{dydt} field contains the derivative vector for the system of equations for the ordinary differential solver. This is created internally from the continuous variables in the processes currently working continuously (i.e., on the continuous event list). The @scheme[get-variable-dt] and @scheme[set-variable-dt!] functions retrieve and store the derivatives of variables to/from the derivative vector when they are controlled by the integration loop. [Note that this field is false, @scheme[#f], when the integration loop is not running.]
The @schemefont{state-changed?} field is true, @scheme[#t], if the system of equations for the ordinary differential solver has changed during a time step and false, @scheme[#f], otherwise. This field is set internally.
The @schemefont{max-step-size} field is the limit of the step size for the ordinary differential equation solver. The default is @scheme[+inf.0], which means no maximum step size. This field may be set by the user.
The following sections show how to use these fields to control the integration loop.
@section{Example---Fixed Step Size---Furnace Model 2a}
If you look at the output of Model 2 above, you will notice a distinct stairstep pattern in the temperature history of the ingots. You will also notice that it overshoots the maximum temperature as much as 20 degrees. This is due to the nature of the near linear temperature function in conjunction with the default control strategy used by the integration loop. The default control strategy will dynamically adjust the step size while keeping the error estimate within an acceptable limit (by default this is @math{10}@superscript{@math{-6}}). Unfortunately, this does work well in this case.
One approach to fixing the model is to use a fixed step size in the integration loop.
In this example, we modify Model 2 to use a fixed step size. In this case, we will use a fixed step size of @math{1} minute.
We set the step size (in the current simulation environment) using the following statement:
@schemeblock[
(current-simulation-step-size (/ 1.0 60.0))
]
where @scheme[(/ 1.0 60.0)] is the step size. [The time unit being used is hours, so @math{1/60} of an hour is one minute.]
Next, we make the integration loop use a fixed step size by setting the @schemefont{control} field (in the current simulation environment) to false, @scheme[#f], as follows:
@schemeblock[
(current-simulation-control #f)
]
These are the only two changes to the model and are included in the updated @schemefont{initialize} function below.
@schemeblock[
(define (initialize)
(code:comment " Set continuous simulation settings")
(current-simulation-step-size (/ 1.0 60.0))
(current-simulation-control #f)
(code:comment " ---")
(set! total-ingots 0)
(set! wait-time (make-variable))
(set! heat-time (make-variable))
(set! leave-temp (make-variable))
(set! pit (make-resource n-pits))
(set! furnace-set (make-set))
(accumulate (variable-history (set-variable-n furnace-set)))
(tally (variable-statistics wait-time))
(tally (variable-statistics heat-time))
(tally (variable-statistics leave-temp))
(tally (variable-history leave-temp))
(schedule #:at 0.0 (scheduler))
(schedule #:at end-time (stop-sim)))
]
The following is the output of this model. Note that the temperature plots are much smoother.
@image["scribblings/images/model-2a-plot-1.gif"]
...
@image["scribblings/images/model-2a-plot-2.gif"]
@verbatim{
Report after 720.0 Simulated Hours - 479 Ingots Processed
-- Ingot Waiting Time Statistics --
Mean Wait Time = 0.3972748535107707
Variance = 0.8502356452674397
Maximum Wait Time = 5.79999999999967
-- Ingot Heating Time Statistics --
Mean Heat Time = 6.877834613068176
Variance = 1.0551156660902024
Maximum Heat Time = 10.813257594378115
Minimum Heat Time = 4.7386143153090075
-- Final Temperature Statistics --
Mean Leave Temp = 901.3208708693385
Variance = 3389.992335826624
Maximum Leave Temp = 1000.3203098942319
Minimum Leave Temp = 801.994264590762}
@image["scribblings/images/model-2a-plot-3.gif"]
@verbatim{
-- Furnace Utilization Statistics --
Mean No. of Ingots = 4.584850935267482
Variance = 3.382482163082148
Maximum No. of Ingots = 7
Minimum No. of Ingots = 0}
@image["scribblings/images/model-2a-plot-4.gif"]
@section{Example---Limiting Step Size---Furnace Model 2b}
Another alternative is to allow the ordinary equation solver control function modify the step size, but limit the step size to a specified maximum.
In this example, we modify Model 2 to use a limited step size. In this case, we will use a maximum step size of @math{1} minute.
We set the maximum step size (in the current simulation environment) using the following statement:
@schemeblock[
(current-simulation-max-step-size (/ 1.0 60.0))
]
where @scheme[(/ 1.0 60.0)] is the step size. [The time unit being used is hours, so @math{1/60} of an hour is one minute.]
These is the only change to the model and are included in the updated @schemefont{initialize} function below.
@schemeblock[
(define (initialize)
(code:comment " Set continuous simulation settings")
(current-simulation-max-step-size (/ 1.0 60.0))
(code:comment " ---")
(set! total-ingots 0)
(set! wait-time (make-variable))
(set! heat-time (make-variable))
(set! leave-temp (make-variable))
(set! pit (make-resource n-pits))
(set! furnace-set (make-set))
(accumulate (variable-history (set-variable-n furnace-set)))
(tally (variable-statistics wait-time))
(tally (variable-statistics heat-time))
(tally (variable-statistics leave-temp))
(tally (variable-history leave-temp))
(schedule #:at 0.0 (scheduler))
(schedule #:at end-time (stop-sim)))
]
The following is the output from the model. The temperature plots are omitted since they are nearly identical to Model 2a.
@verbatim{
Report after 720.0 Simulated Hours - 479 Ingots Processed
-- Ingot Waiting Time Statistics --
Mean Wait Time = 0.39728613689438286
Variance = 0.8495464450741739
Maximum Wait Time = 5.8033836666662495
-- Ingot Heating Time Statistics --
Mean Heat Time = 6.8777413568066175
Variance = 1.0551344984370417
Maximum Heat Time = 10.816121927711492
Minimum Heat Time = 4.7386143153090075
-- Final Temperature Statistics --
Mean Leave Temp = 901.3138280211458
Variance = 3391.192401219392
Maximum Leave Temp = 1000.4429841574237
Minimum Leave Temp = 800.9274096651966}
@image["scribblings/images/model-2b-plot-1.gif"]
@verbatim{
-- Furnace Utilization Statistics --
Mean No. of Ingots = 4.584788893949024
Variance = 3.382567115580496
Maximum No. of Ingots = 7
Minimum No. of Ingots = 0}
@image["scribblings/images/model-2b-plot-2.gif"]
@section{Example---Furnace Model 3}
Furnace Model 3 extends the furnace model (Model 2b) by including a continuous model of the furnace itself. This shows how continuous models may use cooperating continuous processes for more complex simulations.
@schememod[
racket/base
(code:comment " Model 3 - Continuous Simulation Model")
(require (planet williams/simulation/simulation-with-graphics))
(code:comment "Simulation Parameters")
(define end-time 720.0)
(define n-pits 7)
(define initial-furnace-temp 1000.0)
(code:comment "Data collection variables")
(define total-ingots 0)
(define wait-time #f)
(define heat-time #f)
(define leave-temp #f)
(code:comment "Model Definition")
(define random-sources (make-random-source-vector 4))
(define furnace-set #f)
(define furnace-temp #f)
(define pit #f)
(define (scheduler)
(for ((i (in-naturals)))
(schedule #:now (ingot i))
(wait (random-exponential (vector-ref random-sources 0) 1.5))))
(define-process (furnace)
(set! furnace-temp (make-continuous-variable initial-furnace-temp))
(work/continuously
(set-variable-dt!
furnace-temp (* (- 2500.0 (variable-value furnace-temp)) 0.05))))
(define-process (ingot i)
(let* ((initial-temp (random-flat (vector-ref random-sources 1) 100.0 200.0))
(heat-coeff (+ (random-gaussian
(vector-ref random-sources 2) 0.05 0.01) 0.07))
(final-temp (random-flat (vector-ref random-sources 3) 800.0 1000.0))
(current-temp (make-continuous-variable initial-temp))
(arrive-time (current-simulation-time))
(start-time #f))
(with-resource (pit)
(set-variable-value! wait-time (- (current-simulation-time) arrive-time))
(set-insert! furnace-set self)
(set! start-time (current-simulation-time))
(work/continuously
until (>= (variable-value current-temp) final-temp)
(set-variable-dt!
current-temp
(* (- (variable-value furnace-temp) (variable-value current-temp))
heat-coeff)))
(set-variable-value! heat-time (- (current-simulation-time) start-time))
(set-variable-value! leave-temp (variable-value current-temp))
(set-remove! furnace-set self))
(when (variable-history current-temp)
(write-special (history-plot (variable-history current-temp)
(format "Ingot ~a Temperature History" i)))
(newline))
(set! total-ingots (+ total-ingots 1))))
(define (stop-sim)
(printf "Report after ~a Simulated Hours - ~a Ingots Processed~n"
(current-simulation-time) total-ingots)
(printf "~n-- Ingot Waiting Time Statistics --~n")
(printf "Mean Wait Time = ~a~n" (variable-mean wait-time))
(printf "Variance = ~a~n" (variable-variance wait-time))
(printf "Maximum Wait Time = ~a~n" (variable-maximum wait-time))
(printf "~n-- Ingot Heating Time Statistics --~n")
(printf "Mean Heat Time = ~a~n" (variable-mean heat-time))
(printf "Variance = ~a~n" (variable-variance heat-time))
(printf "Maximum Heat Time = ~a~n" (variable-maximum heat-time))
(printf "Minimum Heat Time = ~a~n" (variable-minimum heat-time))
(printf "~n-- Final Temperature Statistics --~n")
(printf "Mean Leave Temp = ~a~n" (variable-mean leave-temp))
(printf "Variance = ~a~n" (variable-variance leave-temp))
(printf "Maximum Leave Temp = ~a~n" (variable-maximum leave-temp))
(printf "Minimum Leave Temp = ~a~n" (variable-minimum leave-temp))
(write-special (history-plot (variable-history leave-temp)
"Final Temperature Histogram"))
(newline)
(printf "~n-- Furnace Utilization Statistics --~n")
(printf "Mean No. of Ingots = ~a~n"
(variable-mean (set-variable-n furnace-set)))
(printf "Variance = ~a~n"
(variable-variance (set-variable-n furnace-set)))
(printf "Maximum No. of Ingots = ~a~n"
(variable-maximum (set-variable-n furnace-set)))
(printf "Minimum No. of Ingots = ~a~n"
(variable-minimum (set-variable-n furnace-set)))
(write-special (history-plot (variable-history (set-variable-n furnace-set))
"Furnace Utilization History"))
(newline)
(stop-simulation))
(define (initialize)
(current-simulation-max-step-size (/ 1.0 60.0))
(set! total-ingots 0)
(set! wait-time (make-variable))
(set! heat-time (make-variable))
(set! leave-temp (make-variable))
(set! pit (make-resource n-pits))
(set! furnace-set (make-set))
(accumulate (variable-history (set-variable-n furnace-set)))
(tally (variable-statistics wait-time))
(tally (variable-statistics heat-time))
(tally (variable-statistics leave-temp))
(tally (variable-history leave-temp))
(schedule #:now (furnace))
(schedule #:at 0.0 (scheduler))
(schedule #:at end-time (stop-sim)))
(define (run-simulation)
(with-new-simulation-environment
(initialize)
(start-simulation)))
(run-simulation)
]
The following is the output from the model.
@verbatim{
Report after 720.0 Simulated Hours - 480 Ingots Processed
-- Ingot Waiting Time Statistics --
Mean Wait Time = 0.0037962597019466577
Variance = 0.0023396090999534486
Maximum Wait Time = 0.9183398430049294
-- Ingot Heating Time Statistics --
Mean Heat Time = 3.3535272061684376
Variance = 0.43028532294776056
Maximum Heat Time = 8.15505100471844
Minimum Heat Time = 2.340282030808339
-- Final Temperature Statistics --
Mean Leave Temp = 902.1065780875852
Variance = 3386.4605258647352
Maximum Leave Temp = 1002.3956405579635
Minimum Leave Temp = 802.2687561089243}
@image["scribblings/images/model-3-plot-1.gif"]
@verbatim{-- Furnace Utilization Statistics --
Mean No. of Ingots = 2.240260843961867
Variance = 2.223332494318446
Maximum No. of Ingots = 7
Minimum No. of Ingots = 0}
@image["scribblings/images/model-3-plot-2.gif"]