Chapter 10

Continuous Simulation Models

The PLT Scheme 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 defining the continuous model.

The PLT Scheme Simulation Collection uses the ordinary differential equation (ODE) solver from the PLT Scheme Science Collection[13] 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.

10.1 Continuous Variables

Continuous variables are used in continuous simulation models. Each continuous variable is associated with the process instance that creates it and when that process instances executes a work/continuously statement, its continuous variables are added to the state vector for the ODE solver.

Continuous variables are variables and have the same data collection facilities available.

make-continuous-variable

Function:

(make-continuous-variable initial-value)
(make-continuous-variable)

Contract:

(case-> (-> real? variable?)
        (-> variable?))

This function returns a new continuous variable with the given initial-value. If no initial-value is specified, 'uninitialized is used. The continuous variable is associated with the current process. The make-continuous-variable call is invalid outside the context of a process.

10.2 Simulation Control (Continuous)

The PLT Scheme 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 macro. This is analogous to the wait/work function in a discrete model, but also specifies the differential equations and, optionally, the terminating conditions, for the continuous model.

work/continuously

Macro:

(work/continuously
  [until expr]
  body ...)

The work/continuously defines a continuous model within a process instance. The body expressions should compute the values of the derivatives of the relevant continuous variables in the process and set them using the set-variable-dt! function. Note that side effects in the body expressions should be avoided since the integration loop will typical evaluate the body expressions multiple times for each time step and may even regress in time when the error estimates in the ODE control object are too high.

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

10.2.1 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[1].

; Model 2 - Continuous Simulation Model
(require (planet "simulation-with-graphics.ss"
                 ("williams" "simulation.plt")))
(require (planet "random-distributions.ss"
                 ("williams" "science.plt")))

;; 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)
  (let loop ((i 0))
    (schedule now (ingot i))
    (wait (random-exponential (vector-ref random-sources 0) 1.5))
    (loop (+ i 1))))

(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)
      (if (= (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))
    (if (variable-history current-temp)
        (printf
         "~a~n"
         (history-plot (variable-history current-temp)
                       (format "Ingot ~a Temp History" i))))
    (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))
  (printf "~a~n"
          (history-plot (variable-history leave-temp)
                        "Final Temperature Histogram"))
  (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)))
  (printf
   "~a~n"
   (history-plot (variable-history (set-variable-n furnace-set))
                 "Furnace Utilization History"))
  (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 end-time) (stop-sim))
  (schedule (at 0.0) (scheduler)))

(define (run-simulation)
  (with-new-simulation-environment
    (initialize)
    (start-simulation)))

>(run-simulation)

. . .

[continuous-Z-G-2.gif]

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

10.3 Simulation Environment (Continuous)

If you look at the output of Model 2 above you will notice a distinct stair-step pattern in the temperature history of the ingots. You will also notice that it overshoots the maximum temperature by 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 with an acceptable limit (by default this is 10^-6). Unfortunately, this doesn't work well in this case.

10.3.1 Fixed Step Size

10.3.1.1 Example - Furnace Model 2a

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 $1$ minute. Next, we make the integration loop use a fixed step size by setting the ODE control object to \texttt{\#f}. \scm{ (current-simulation-control #f)} These are the only changes to the model and are included in the updated \verbinitialize+ 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 end-time) (stop-sim))
  (schedule (at 0.0) (scheduler)))

The following is the output of this model. Note that the temperature plots are much smoother.

>(run-simulation)

. . .

[continuous-Z-G-6.gif]

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

10.3.2 Limiting Step Size

Another alternative is to allow the ODE control object to modify the step size, but limit it.

10.3.2.1 Example - Furnace Model 2b

In this example we modify Model 2 to use a limited step size. In this case we will limit the step size to a maximul 1 minute.

We set the step size limit (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 maximum step size. [The time unit being used is hours, so $1/60$ of an hour is $1$ minute. This is the only change to the model and is included in the updated \verbinitialize+ 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 end-time) (stop-sim))
  (schedule (at 0.0) (scheduler)))

The following is the output from the model. The temperature plots are omitted since they are nearly identical to Model 2a.

>(run-simulation)

. . .

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

10.4 Furnace Model 3

Model 3 extents model by including a continuous model of the furnace.

; Model 3 - Continuous Simulation Model
(require (planet "simulation-with-graphics.ss"
                 ("williams" "simulation.plt")))
(require (planet "random-distributions.ss"
                 ("williams" "science.plt")))

;; 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)
  (let loop ((i 0))
    (schedule now (ingot i))
    (wait (random-exponential (vector-ref random-sources 0) 1.5))
    (loop (+ i 1))))

(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))
;    (if (= (modulo i 100) 0)
;        (accumulate (variable-history current-temp)))
    (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))
    (if (variable-history current-temp)
        (printf "~a~n"
                (history-plot (variable-history current-temp)
                              (format "Ingot ~a Temperature History" i))))
    (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))
  (printf "~a~n"
          (history-plot (variable-history leave-temp)
                        "Final Temperature Histogram"))
  (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)))
  (printf "~a~n"
          (history-plot (variable-history (set-variable-n furnace-set))
                        "Furnace Utilization History"))
  (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 (at end-time) (stop-sim))
  (schedule (at 0.0) (scheduler))
  (schedule now (furnace)))

(define (run-simulation)
  (with-new-simulation-environment
    (initialize)
    (start-simulation)))

(run-simulation)

>(run-simulation)
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.1065780875851
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