Making Connections

STELLA/iThink version 10 introduces several new array features, including simplified and more powerful Apply-To-All equations that are designed to reduce the need to specify equations for every individual element.

Dimension names are optional

When an equation is written using other array names, the dimension names are not normally needed.  For example, given arrays A, B, and C, each with the dimensions Dim1 and Dim2, A can be set to the sum of B and C with this equation:

B + C

Dimension names are still needed when the dimensions do not match.  For example, to also add in the first 2-dimensional slice of the 3-dimensional array D[Dim1, Dim2, Dim3], the equation becomes:

B + C + D[Dim1, Dim2, 1]

The wildcard * is optional

When an array builtin is used, the * is normally not needed.  For example, to find the sum of the elements of a 2-dimensional array A[Dim1, Dim2] requires this equation:

SUM(A)

If, however, the sum of only the first column of A is desired, the * is still needed:

SUM(A[*, 1])

Simplified array builtins

There are five array builtins:  SIZE, SUM, MEAN, STDDEV, and RANK.  In addition, the MIN and MAX functions have been extended to take either one or two array arguments.  All but RANK can also be applied to queues and conveyors.

SUM, MEAN, and STDDEV all work in a similar way (see examples of SUM above).

Using the MAX function, it is possible to find the maximum value in array A,

MAX(A)

the maximum value in array A, or zero if everything is negative,

MAX(A, 0)

or the maximum across two arrays A and B,

MAX(A, B)

MIN works the same way, but finds the minimum.

The SIZE function requires an array parameter, but within an array, the special name SELF can be used to refer to the array whose equation is being set.  In addition, wildcards can be used to determine the size of any array slice.  In the equation for array A[Dim1, Dim2],

SIZE(SELF)

gives the total number of elements in array A while

SIZE(SELF[*, 1])

gives the size of the first dimension of A, i.e., the number of elements – or rows – in the first column.  Likewise,

SIZE(SELF[1, *])

gives the size of the second dimension of A, i.e., the number of elements – or columns – in the first row.

Since RANK returns the index of the element with the given rank, it can also be used to find the index of the minimum element (using rank 1) or the maximum element (using rank SIZE(array)).  Given array A[Dim1, Dim2], the index of the minimum element in the first row can be found with the equation:

RANK(A[1, *], 1)

However, to find the minimum element in the entire array, use:

RANK(A, 1)

This returns a single index that can be mapped to an array element using the special parentheses subscripting:

A(RANK(A, 1))

will be the value of the minimum element in A, i.e, the same value as MIN(A).  However, if array B has the same dimensions as A (i.e., for this example, B[Dim1, Dim2]), the value of the element in B that corresponds to the minimum element in A is found with:

B(RANK(A, 1))

Accessing elements of queues and conveyors

Use an array subscript to access an element of a queue or conveyor.  The indices start on the outflow side (at 1) and increase toward the inflow side (up to SIZE(queue) or SIZE(conveyor)).  This allows the entire contents of a queue or conveyor to be assigned to an array allowing additional calculations, for example, a weighted average.  Given a conveyor named Lag, a new array weighted_by_time[Slat] can be created with the equation:

(Slat*DT)*Lag[Slat]

Note the subscript is required for the conveyor.  Otherwise, the total value of the conveyor will be used.  Note also that the size of the dimension Slat must be at least large enough to hold all of the conveyor elements (the remaining elements in weighted will be set to zero).  The value of Slat*DT is the amount of time remaining before the material in that slat exits the conveyor.

A converter, average_latency, which is the average time remaining for the contents to exit (a weighted mean), can now be defined with the equation:

SUM(weighted_by_time)/Lag

Transposition

It is sometimes helpful to transpose an array.  To facilitate this, the ' (apostrophe) operator was added.  Given arrays A[Dim1, Dim2, Dim3] and B[Dim3, Dim2, Dim1], the array A can be set equal to B transposed with the following equation:

B '

Note that a space is required between the array name and the apostrophe.  This is equivalent to the following equation that uses dimension names:

B[Dim3, Dim2, Dim1]

This is especially helpful for square matrices or other arrays that use the same dimension name many times.  Given arrays C[Dim, Dim, Dim] and D[Dim, Dim, Dim], the array C can be set equal to D transposed with the following equation, which reverses all the dimensions:

D '

This is equivalent to the following equation that uses the new positional dimension names:

D[@3, @2, @1]

Within a subscript, the @ operator can be followed by an integer that represents the dimension position in the array whose equation is being set.  In the example above, @3 represents the third dimension name of A.  This is particularly useful if straight transposition is not needed and all the dimension names are the same.  For example,

D[@2, @1, @3]

flips the first two dimensions of D (when assigning to A) while leaving the third alone.

Subscript expressions

Subscripts can contain any valid expression.  Given an array A and a variable x, an element at a variable index that is one more than twice x can be accessed with:

A[2*x + 1]

Element labels can also appear within these expressions.

In Apply-To-All arrays, dimension names can be used.  The following equation sets the values in array A[Dim1] to every even-indexed elements in array B[Dim1], filling the second half of A with zeroes:

B[2*Dim1]

Dimension names can also be used outside subscripts.  The following equation slides the elements of B up one position in A, placing 10 in the first element of A (without the IF, the first element would contain 0).

IF Dim1 = 1 THEN 10 ELSE B[Dim1 - 1]

Even if Dim1 is labeled, it must be compared to the numeric index 1 in the IF expression because element labels can only be used within a subscript.  Note that numeric indices are always valid for any array dimension, even if it is labeled.

Array Ranges

A range of an array can be specified using the range operator : (colon), which takes a lower bound on the left and an upper bound on the right (e.g., 1:10 means “from 1 to 10”).  Just as wildcards allow control over which dimensions to include, ranges control which range of elements to include in each dimension.  For example, the follow equation sums the top-left 3×4 rectangle of array A[Dim1, Dim2]:

SUM(A[1:3, 1:4))

We hope you find these new array capabilities useful in your modeling work and welcome any comments and suggestions.

There are times when it is helpful to calibrate, or fit, your model to historical data. This capability is not built into the iThink/STELLA program, but it is possible to interface to external programs to accomplish this task. One generally available program to calibrate models is PEST, available freely from www.pesthomepage.org. In this blog post, I will demonstrate how to calibrate a simple STELLA model using PEST on Windows. Note that this method relies on the Windows command line interface added in version 9.1.2 and will not work on the Macintosh. The export to comma-separated value (CSV) file feature, added in version 9.1.2, is also used.

The model and all files associated with its calibration are available by clicking here.

The Model

The model being used is the simple SIR model first presented in my blog post Limits to Growth. The model is shown again below. There are two parameters: infection rate and recovery rate. Technically, the initial value for the Susceptible stock is also a parameter. However, since this is a conserved system, we can make an excellent guess as to its value and do not need to calibrate it.

The Data Set

We will calibrate this model to two data sets. The first is the number of weekly deaths caused by the Hong Kong flu in New York City over the winter of 1968-1969 (below).

The second is the number of weekly deaths per thousand people in the UK due to the Spanish flu (H1N1) in the winter of 1918-1919 (shown later).

In both cases, I am using the number of deaths as a proxy for the number of people infected, which we do not know. This is reasonable because the number of deaths is directly proportional to the number of infected individuals. If we knew the constant of proportionality, we could multiply the deaths by this constant to get the number of people infected.

Read more…

After reading Karim Chichakly’s recent post on Integration Methods and DT, I was reminded that delta time (DT) has always been a tricky modeling concept for me to grasp.   Beginning modelers don’t usually need to think about changing DT since STELLA and iThink set it to a useful default value of 0.25.   But once you progress with your modeling skills, you might consider the advantages and risks of playing with DT.

The DT setting is found in the Run Specs menu.

By definition, system dynamics models run over time and DT controls how frequently calculations are applied each unit of time.  Think of it this way, if your model was a movie, then DT would indicate the time interval between still frames in the strip of movie film.  For a simulation over a period of 12 hours, a DT of 1/4 (0.25) would give you a single frame every 15 minutes.  Lowering the DT to 1/60 would give a frame every minute.   The smaller the DT is, the higher the calculation frequency (1/DT).

Beware of the Extremes

A common tendency for modelers is to set the calculation frequency too high.  Without really thinking too hard about it, more data seems to imply a higher quality model – just like more frames in movie film make for smoother motion.  If your model calculates more data for every time unit, its behavior will begin to resemble the behavior of a smoothly continuous system.  But a higher frequency of calculations can greatly slow down your model’s run performance and more data does not directly translate to a better simulation.

Beware of Discrete Event Models

Another situation where DT can often lead to unexpected behavior is with models that depend on discrete events.   My eyes were opened to this when I attended one of isee’s workshops taught by Corey Peck and Steve Peterson of Lexidyne LLC.

One of the workshop exercises involved a simple model where the DT is set to the default 0.25, the inflow is set to a constant 10, and the outflow is set to flush out the stock’s contents as soon as it reaches 50.   This is how the model’s structure and equations looked:

Stock = 0

inflow = 10

outflow = IF Stock >= 50 THEN 50 ELSE 0

I would have expected the value of the stock to plunge to zero after it reached or exceeded 50, but this graph shows the resulting odd saw-tooth pattern.

The model ends up behaving like a skipping scratched record, in a perpetual state of never progressing far enough to reach the goal of zero.  (Click here to download the model.)

What is happening in the model?  In the first DT after the stock’s value reaches exactly 50, the outflow sets itself to 50 in order to remove the contents from the stock. So far so good, but now the DT gotcha begins to occur.   Since the outflow works over time, its value is always per time.  To get the quantity of material that actually flowed, you must multiply the outflow value (or rate) by how long the material was flowing.  When DT is set to 0.25,  the material flows 0.25 time units each DT.  Hence, the quantity of material removed from the stock is 50*0.25 = 12.50.

Suddenly we are in a situation where only 12.50 has been removed from the stock but the stock’s value is now less than 50.  Since the stock is no longer greater than or equal to 50, the outflow sets itself back to 0 and never actually flushes out the full contents of the stock. 

So what do we do?  One solution to this problem would be to use the PULSE built-in to remove the full value from the stock.   Here’s what the equation for the outflow would look like:

outflow = IF Stock >= 50 THEN PULSE(Stock) ELSE 0

(Note: This option will only work using Euler’s integration method.)

Further Reading

STELLA and iThink have great help documentation on DT.  The general introduction provides a good explanation of how DT works. The more advanced DT Situations Requiring Special Care section focuses more on artifactual delays and the discrete model issues mentioned in this post.  Delta time and resulting model behaviors are reminders that system dynamics models run over time, but they achieve this by applying numerous discrete calculations in order to simulate the smooth behavior of actual systems.

The simulation engine underlying STELLA® and iThink® uses numerical integration.  Numerical integration differs from the integration you may have learned in Calculus in that it uses algorithms that approximate the solution to the integration.  The two approximations currently available are known as Euler’s method and the Runge-Kutta method.  All algorithms require a finite value for DT, the integration step-size, rather than the infinitesimally small value used in Calculus.  On the surface, it may seem that the smaller DT is, the more accurate the results, but this turns out not to be true.

Compound Interest:  Euler’s Method over Runge-Kutta

To introduce Euler’s method, let’s take a look at the simple problem of compound interest.  If we have $100 that we invest at 10% (or 0.1) compounded annually, we can calculate the interest after N years by adding in the interest each year and recalculating:

1st year:  interest = $1000 × 0.1 = $100; Balance = 1000 + 100 = $1100
2nd year: interest = $1100 × 0.1 = $110; Balance = 1100 + 110 = $1210
3rd year:  interest = $1210 × 0.1 = $121; Balance = 1210 + 121 = $1331

And so on up to year N.  We have just seen the essence of how Euler’s method works.  It calculates the new change in the stock for this DT (in this case, interest) and then adds that to the previous value of the stock (Balance) to get the new value of the stock.  In this example, DT = 1 year.

By noticing we always add the existing balance in, we can instead just multiply the previous year’s balance by 1 + rate = 1 + 0.1 = 1.1:

1st year:  Balance = $1000 × 1.1 = $1100
2nd year: Balance = $1100 × 1.1 = $1210
3rd year:  Balance = $1210 × 1.1 = $1331

And so on up to year N. We can further generalize by noticing we are multiplying by 1.1 N times and thus arrive at the compound interest formula:

Balance = Initial_Balance*(1 + rate)^N

Checking this, we find our Balance at the end of year 3 is 1000*1.1^3 = $1331.  In the general case of the formula, rate is the fractional interest rate per compounding period and N is the number of compounding periods (an integer).  In our example, the compounding period is one year, so rate is the annual fractional interest rate and N is the number of years.  However, if interest is compounded quarterly (four times a year), the interest rate has to be adjusted to a per quarter rate by dividing by 4 (so rate = 0.1/4 = 0.025) and N must be expressed as the number of quarters (N = number of years*4 = 3*4 = 12 for the end of year 3).  We can use this formula in our model to test the accuracy of Euler’s method.  Note that for quarterly compounding, we would set DT = 1/4 = 0.25 years.

To explore the differences between Euler’s and Runge-Kutta, the following structure will be used for all of the examples in this post.  This structure models the compound interest problem outlined above.

The equations change for each example and can be seen in the individual model files (accessed by clicking here).  For this example, the actual value is calculated using the compound interest formula, Initial_Balance*(1 + rate)^TIME.  The approximated value is calculated by integrating rate*Approx_Balance (into Approx_Balance).

In addition to the actual and approximate values, three errors are also calculated across the model run:  the maximum absolute error, the maximum relative error, and the root-mean-squared error (RMSE).  The absolute error is:

ABS(Actual_Balance – Approx_Balance)

The relative error is:

absolute_error/ABS(Actual_Balance)

and is usually expressed as a percentage.  The RMSE is found by averaging the values of the absolute error squared, and then taking the square root of that average.

Read more…

Conveyors are useful model elements for representing pipelines or processes that take a certain amount of time to complete.  However, adding a leakage flow to a conveyor can make it difficult to initialize a model in steady-state.  The following discussion will explain how to initialize conveyors with leakage in steady-state.  Please refer to the model structure below while reading this discussion.

These additional variables will be also used:

transit_time = TRANSTIME(conveyor)
conveyor_length = transit_time/DT
leakage_fraction = the user-specified leakage fraction

Linear Leakage

The default leakage is linear in behavior.  The total amount that leaks across the length of the conveyor is directly proportional to the inflowing amount.  The leakage fraction is the constant of proportionality.  Thus, the fraction of inflowing material that makes it to the conveyor’s outflow is exactly

1 – leakage_fraction

Given the sample model structure above, to achieve equilibrium, conveyor_outflow must equal outflow.  For this to happen, we need to set the inflow as follows:

inflow = outflow/(1 – leakage_fraction)

The conveyor’s steady-state value is then:

conveyor = transit_time*inflow – (conveyor_length – 1)*leakage*DT/2

where the initial value of leakage is:

leakage = leakage_fraction*inflow

This must be calculated outside the program and entered as a constant into the conveyor as conveyors cannot be given equations (they can, however, be set to a the value of a single converter, but you must be careful how you calculate this to avoid circularity).

Exponential Leakage

Optionally, leakage can be made exponential.  The amount that leaks each DT is proportional to the amount remaining in the conveyor.  In this case, the leakage fraction is the fraction that leaks each unit of time so, for long conveyors, a lot of material can leak away.  Given the transit time, the fraction of inflowing material that makes it to the conveyor’s outflow is approximately

1 – (1 – leakage_fraction)^transit_time

Given the sample model structure above, to achieve equilibrium, conveyor_outflow must equal outflow.  For this to happen, we need to set the inflow as follows:

per_dt_no_leak = 1 – leakage_fraction*DT
inflow = outflow/(per_dt_no_leak^conveyor_length)

For steady-state, the conveyor itself must then be set as follows:

conveyor = (inflow*DT)*(1 – per_dt_no_leak^conveyor_length)/(1 – per_dt_no_leak)

I generally do not use modules to build very small models (only a couple of stocks and flows), which may then lead me to use sectors as the model grows because they are very convenient.  By the time I have three sectors, though, it starts to become clear that I should have used modules.  I will then need to convert my sector-based model into a module-based model.  Historically, I also have a number of sector-based models that are crying to be module-based.

Converting from sectors to modules is not very difficult:

1. Make sure there are no connections or flows between sectors.  Replace any of these with ghosts in the target sector.
2. In a new model, create one module for every sector.
3. Copy and paste the structure from each sector into its corresponding module.
4. Connect the modules:  At this point, the model structure has been rearranged into modules, but none of the modules are connected.  The ghosts that were in the sectors became real entities when they were pasted into the modules.  Go back to identify all of these connections and reconnect them in the module-based model.

Stepping Through a Sample Model

Let’s walk through an example.  A small sector-based model is shown below (and is available by clicking here).

This model violates what I would call good sector etiquette:  there are connectors that run between the sectors.  This is often useful in a small model such as this because it makes the feedback loops visible.  However, in a larger model, this can lead to problems such as crossed connections and difficulty in maintaining the model because sectors cannot be easily moved.

Read more…

This is the last installment of a four-part series.  The first three parts can be accessed by clicking on the links below.
Methods for Using Arrays Effectively
Modeling a Watershed with Arrays
Modeling Customers Switching Between Brands

Generalizing the Model

When I showed Steve Peterson (at Lexidyne) my brand switching model, he told me there is a more general version that separates the customer loss fraction from the fraction won by another competitor.  This has been presented in Pharmaceutical Product Strategy by Mark Paich, Corey Peck, and Jason Valant.

In my original formulation, the switching probability matrix was the product of these two variables.  However, in many practical cases, the data available comes from two different places and reflects these two separate components.  The revised model structure is shown below.

Instead of one composite switching probability, this model uses a switching out probability that is distinct from the switching in probability.  The switching out probability is a one-dimensional array that, for each product, contains the fraction of customers lost to rivals every time unit (in our case, month).  A sample for the five brands A, B, C, D, and E appears below.

switching out probability

We can see from this table that Brand B is losing 17% of its customers to rivals each and every month!  Whoever is managing that product had better do something quickly.

The other side of the story has to do with which brand the customers are switching to.  The switching in probability matrix contains, for each brand, the fraction of lost customers that migrate to a rival brand.  Thus, each row of this matrix must add up to one (100% of lost customers).  A sample appears below.

switching in probability

Note the diagonal will always be zero.

We can determine a lot of things from this table.  For example, brand B offers no competition to brand D, brand D is the biggest rival of all the other brands, and brand C is brand D’s biggest rival.

Read more…

This is an update to post published on August 31, 2009.  The attached model was updated to find negative means and an alternate method was included at the end.

I am frequently asked which built-in function gives either the running mean or running standard deviation of a model variable.  Unfortunately, there is no such built-in at this time (no, that is not what MEAN() does).

Luckily, however, we can replicate the behavior we desire from built-in functions by creating a reusable module.  I can create a module that calculates a running average and a running standard deviation from any model variable.

When building a reusable module component, it is important to carefully define what the input to the module will be (i.e., what are the parameters to the built-in function) and what the output of the module will be (i.e., what is the result or return value of the built-in function).  In this particular case, the input will be the variable whose running average or running standard deviation we wish to find.  There are two outputs:  the running average and the running standard deviation.  Note we do not have to use both outputs all the time.

Thus, our new module can be used as shown below:

Note the name of the module was chosen to give a meaningful context to the running mean and standard deviation variables, which have fixed names defined inside the reusable module.  As in this example, it is always a good idea to give the module outputs general names that make sense when qualified by a context (the module name).

The reusable module itself was built and tested in iThink, and can also be used in STELLA.  The input parameter was given an equation to allow the model to be completely tested and debugged before being reused.  The model appears below and can be downloaded by clicking here.

Note the input to the module is named value.  After importing the module, this will need to be assigned to the variable in question, Cash in the above example.  This can be done from outside the module by right-clicking on Cash and choosing “Module->Assign to”, or right-clicking on value and choosing “Module->Assign Input to”.  The outputs can be assigned in a similar way, or the Ghost tool can be used.

This method, while relatively easy to understand, does accurately compute the standard deviation when the mean of the running sum of squares is close in magnitude to the running mean squared.  An alternate method that does not suffer this problem was developed by Welford in 1962 and is implemented in the model that can be downloaded by clicking here.

Finally, I am including a simple reusable module that finds the maximum value of a model variable across the entire run of a simulation.  It can be downloaded by clicking here.  It uses a stock to hold the maximum value seen so far, and takes advantage of the fact that uniflows cannot be negative.  It is used the same way as the running mean and standard deviation module, but only has one output called maximum.

For more information about modules, consult the iThink and STELLA help files.  These on-line resources are also available:

Using Modules Webinar

Module FAQs