11 Continuous Simulation Models

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 wait/continuously statement that specifies the differential equations for the continuous model and, optionally, the termination condition for the model.

11.1 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 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.

The variable field continuous? specifies whether or not a variable is a continuous variable.

The variable field state-index is only used for continuous variables. When the process instance owning the continuous variable is in the PROCESS-WORKING-CONTINUOUSLY state (that is, is currently in a 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 get-variable-value) are retrieved from the state vector. The state-index field is -1 for discrete variables or when a continuous variable is not being controlled by the integration loop.

11.2 Simulation Control (Continuous)

The simulation collection main loop (i.e., the start-simulation function) implements the integration loop internally.

A continuous simulation model is implemented within a process using the work/continuously statement. This is analogous to the wait/work statement in a discrete simulation model, but specifies the differential equations and, optionally, the terminating condition for the continuous model.

The expr in the until clause is evaluated at the start of each time step. If the expr is true, the process drops out of the work/continuously and resumes execution as a discrete event simulation model.

11.3 Example—Furnace Model 2

This combined discrete-continuous simulation model is derived from an example in Introduction to Combined Discrete-Continuous Simulation Using SIMSCRIPT II.5 by Abdel-Moaty M Fayek [Fayek02].

#lang racket/base

; Model 2 - Continuous Simulation Model

(require (planet williams/simulation/simulation-with-graphics))

; Simulation Parameters

(define end-time 720.0)

(define n-pits 7)

; Data collection variables

(define total-ingots 0)

(define wait-time #f)

(define heat-time #f)

(define leave-temp #f)

; 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.

...

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

-- Furnace Utilization Statistics --

Mean No. of Ingots    = 4.687638930905194

Variance              = 3.321384522948101

Maximum No. of Ingots = 7

Minimum No. of Ingots = 0

11.4 Simulation Environment (Continuous)

There are several field in each simulation environment that are specific to continuous models.

The 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 work/continuously statement. These (should) implement the continuous models by computing the appropriate variable derivatives. This field is set internally.

The 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 control field contains the ordinary differential equation control object of false, #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 10-6 and relative error of 0.0 with respect to the solution. This field may be set by the user—in particular, set this field to false, #f, for a fixed step size.

The 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 step field contains the ordinary equation stepper function. The stepper function advances the solution from time t to time t + h for a fixed step size h and estimates the resulting local error. This field is set internally.

The 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 step-size field contains the current (last) step size for the ordinary differential equation solver. If the control field is false, "#f", this specifies the fixed step size; otherwise, it is controlled by the control object. The default step size is 10-6 (only applicable when the control field is false, #f). This field may be set by the user.

The dimension field contains the dimensionality of the system of equations for the ordinary differential solver. This is set internally.

The 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 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, #f, when the integration loop is not running.]

The 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 get-variable-dt and 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, #f, when the integration loop is not running.]

The state-changed? field is true, #t, if the system of equations for the ordinary differential solver has changed during a time step and false, #f, otherwise. This field is set internally.

The max-step-size field is the limit of the step size for the ordinary differential equation solver. The default is +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.

11.5 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 10-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 1 minute.

We set the step size (in the current simulation environment) using the following statement:

(current-simulation-step-size (/ 1.0 60.0))

where (/ 1.0 60.0) is the step size. [The time unit being used is hours, so 1/60 of an hour is one minute.]

Next, we make the integration loop use a fixed step size by setting the control field (in the current simulation environment) to false, #f, as follows:

(current-simulation-control #f)

These are the only two changes to the model and are included in the updated initialize function below.

(define (initialize)

;  Set continuous simulation settings

(current-simulation-step-size (/ 1.0 60.0))

(current-simulation-control #f)

;  ---

(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.

...

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

-- Furnace Utilization Statistics --

Mean No. of Ingots    = 4.584850935267482

Variance              = 3.382482163082148

Maximum No. of Ingots = 7

Minimum No. of Ingots = 0

11.6 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 1 minute.

We set the maximum step size (in the current simulation environment) using the following statement:

(current-simulation-max-step-size (/ 1.0 60.0))

where (/ 1.0 60.0) is the step size. [The time unit being used is hours, so 1/60 of an hour is one minute.]

These is the only change to the model and are included in the updated initialize function below.

(define (initialize)

;  Set continuous simulation settings

(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 #: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.

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

-- Furnace Utilization Statistics --

Mean No. of Ingots    = 4.584788893949024

Variance              = 3.382567115580496

Maximum No. of Ingots = 7

Minimum No. of Ingots = 0

11.7 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.

#lang racket/base

;  Model 3 - Continuous Simulation Model

(require (planet williams/simulation/simulation-with-graphics))

; Simulation Parameters

(define end-time 720.0)

(define n-pits 7)

(define initial-furnace-temp 1000.0)

; Data collection variables

(define total-ingots 0)

(define wait-time #f)

(define heat-time #f)

(define leave-temp #f)

; 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.

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

-- Furnace Utilization Statistics --

Mean No. of Ingots    = 2.240260843961867

Variance              = 2.223332494318446

Maximum No. of Ingots = 7

Minimum No. of Ingots = 0