11 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 implementing the continuous model.
The PLT Scheme Simulation Collection uses the ordinary differential equation (ODE) solver from the PLT Scheme 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.
| (make-continuous-variable initial-value) → variable? |
| initial-value : (or/c real? (symbols uninitialized)) |
Returns a new continuous variable with the specified initial-value. If no initial-value is specified, 'uninitialized is used. The continuous variable is associated with the current process instance. A make-continuous-variable call is invalid outside of a process.
The variable field continuous? specifies whether 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 and when a continuous variable is not being controlled by the integration loop.
| (variable-dt! variable) → real? |
| variable : variable? |
Returns the value of the derivative of 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 PROCESS-WORKING-CONTINUOUSLY state).
| (set-variable-dt! variable value) → any |
| variable : variable? |
| value : real? |
Sets the value of the derivative of variable to 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 PROCESS-WORKING-CONTINUOUSLY state). This should be used in the work/continuously to implement the differential equations for the continuous model.
11.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 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.
| |||
|
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 typically evaluate the 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 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 scheme/base |
| ; Model 2 - Continuous Simulation Model |
| (require (planet williams/simulation/simulation-with-graphics)) |
| (require (planet williams/science/random-distributions)) |
| ; 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) |
| (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 end-time) (stop-sim)) |
| (schedule (at 0.0) (scheduler))) |
| (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 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.
...
| 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 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.
| 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 scheme/base |
| ; Model 3 - Continuous Simulation Model |
| (require (planet williams/simulation/simulation-with-graphics)) |
| (require (planet williams/science/random-distributions)) |
| ; 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)) |
| (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 (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) |
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 |
