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.

Sarah Boyar accepts Scholarship Award from Joanne Egner

The Barry Richmond Scholarship Award was established in 2007 by isee systems to honor and continue the legacy of its founder, Barry Richmond.  Barry was devoted to helping others become better systems citizens.  Systems citizens are members of a global community that strive to understand the complexities of today’s world and have the informed capacity to make a positive difference.  It was Barry’s mission to make systems thinking and system dynamics accessible to people of all ages, and in all fields.  The award is presented annually at the System Dynamics Society Conference to an individual whose work demonstrates a desire to expand the field of systems thinking or to apply it to current social issues.

Through most of his career, Barry focused on education as the key to spreading systems thinking.  As a teacher and a mentor he dedicated much of his time to developing tools and methodologies for teaching systems thinking.  With this in mind, it was a great pleasure to present this year’s award to Sarah Boyar, a recent graduate of the Masters Program in System Dynamics at Worcester Polytechnic Institute (WPI).

Sarah Boyar and Karim Chichakly enjoy the conference banquet

Sarah presented a portfolio of her work to the scholarship committee.  In particular, an essay about her teaching philosophy resonated with us.  Sarah wrote this piece while taking a seminar in college teaching in order to fulfill her Real World Dynamics course requirement at WPI.  Since she already had plenty of experience as a consultant applying system dynamics to real world situations, Sarah managed to convince the WPI powers-that-be that an essential real world manifestation of system dynamics is the way that it is taught.  This is something Barry would have encouraged and been excited about.

Her essay titled Beliefs About Teaching and Learning begins as follows:

I teach System Dynamics. While I want my students to have some knowledge of system dynamics, most of all I want them to be excited and stimulated by it. I also want them to find it beautiful: I want to teach in such a way that my students find some aspect of beauty in the work, whether it’s through the visual arcs in the model interface, or the precision of algebra in the way we write statements, or the way that system dynamics can ameliorate a social ill that concerns them. I want my students to somehow feel a sense of peace and beauty derived from some aspect of the knowledge I am teaching.

Among Sarah’s aspirations is the desire to teach system dynamics to professionals in other fields, namely lawyers (potential judges) and medical doctors.  Incorporating a systems perspective within both the judicial system and in healthcare could certainly make a positive difference for us all.  Good luck and congratulations Sarah!

Essentially, data can be moved between a database and iThink/STELLA by using a CSV file as a bridge. CSV files are a widely supported file standard for storing table data, and both iThink/STELLA and many database programs are able to read and write to them.

The process of connecting to a database using CSV files as an intermediary

The process can be automated when you use iThink/STELLA’s ability to run models automatically from the command line (Windows only). Most database applications also have command line interfaces, allowing you to create a single macro script that moves data between your model and a database in a single process.

In this post I will use a simple example to demonstrate how to import data from a Microsoft SQL Server database into an iThink model on Windows. The model and all files associated with the import process are available by clicking here. If you don’t have access to Microsoft SQL Server, you can download a free developer’s version called SQL Server Express from the Microsoft web site.

The Model

The model used in this example is a variation of the Beer Game model. The structure shown below represents the ordering process for a simple retailer supply chain.

The model has been set up to import the initial values for On Order with Wholesaler and Unfilled Orders stocks, target inventory and actual customer orders (a graphical function with 21 weeks of data). The source of the imported data is the file named import.csv in the example files.

To set up this example, I manually created the CSV file using the initial model parameters. (Later in this post, you’ll see that this file will be automatically created by the database.) The model has been initialized in a steady state with actual customer orders at a constant level of 4 cases per week over the 21 week period.

Creating the SQL Server Database

Since the model is set up to import data from the import.csv file, we are ready to populate that same CSV file with data from the SQL Server database. But, first we need to create the database and the associated data table.

Using SQL Server Management Studio, create a new database named Simulation.

Dialog in SQL Server when adding new "Simulation" database

Once the database has been created, you can run the full CreateTable.sql script to create a table that will store the data. To open the script, select Open -> File from the File menu and navigate to the CreateTables.sql file. The script is set up to create a table named Import_Data with columns for each variable that is imported into the model. The script will also populate the table with data. It is important to note that the first row of the table contains the variable names exactly as they appear in the retailer supply chain model. Below is a sampling of the commands contained in the script:

CREATE TABLE [dbo].[Import_Data](
[RowID] [int] IDENTITY(1,1) NOT NULL,
[OnOrder] [nvarchar](50) COLLATE SQL_Latin1_General_CP1_CI_AS NULL,
[UnfilledOrders] [nvarchar](50) COLLATE SQL_Latin1_General_CP1_CI_AS NULL,
[TargetInventory] [nvarchar](50) COLLATE SQL_Latin1_General_CP1_CI_AS NULL,
[ActualOrders] [nvarchar](50) COLLATE SQL_Latin1_General_CP1_CI_AS NULL
) ON [PRIMARY]

INSERT INTO Import_Data
(OnOrder,UnfilledOrders,TargetInventory,ActualOrders)
values
(‘On Order with Wholesaler’,'Unfilled Orders’,'target inventory’,'actual customer orders’)
INSERT INTO Import_Data
(OnOrder,UnfilledOrders,TargetInventory,ActualOrders)
values
(8,4,20,4)
INSERT INTO Import_Data
(OnOrder,UnfilledOrders,TargetInventory,ActualOrders)
values
(”,”,”,4)
INSERT INTO Import_Data
(OnOrder,UnfilledOrders,TargetInventory,ActualOrders)
values
(”,”,”,10)

To run the script, select Execute from the Query menu. Note that the data shows an increase in actual customer orders in week 2 from 4 to 10 cases. In the model, actual customer orders is a graphical function which is reflecting order data over time.

Results tab after running CreateTable.sql script

Linking the Database to the Model

After the database has been prepared, the Transfer_Data.bat Windows batch file can be used to link the database to the model. Transfer_Data.bat is a script that pulls data from the database into the CSV file, opens iThink, imports the CSV data and runs the model. The batch file takes advantage of the utility Microsoft provides with SQL Server called Bulk Copy Program (BCP) and iThink’s ability to run models from the command line. You’ll need to make a few edits to the batch file shown below before you can run it.

BCP Simulation.dbo.Import_Data out import.csv -c -t, -U <Login> -P <Password> -S <Server name>

“%PROGRAMFILES%\isee systems\iThink 9.1.4\iThink.exe” -i -r -nq Supply_Chain.itm

Using a text editor such as WordPad or NotePad, replace the placeholders for <Login>, <Password> and <Server name> in the first line to match the authentication for the Simulation database that was previously created. If your SQL Server database is on your local machine, you can delete the –U, -P and –S settings altogether. This command will make use of the BCP utility to copy the contents of the Import_Data table to the import.csv file. Depending on the version of iThink or STELLA you are running (version 9.1.2 or later is required) and its location, you may also need to edit the second line of the batch file.

After the Transfer_Data.bat file has been edited and saved, double-click on the file to run it. Once the data has been retrieved from the database, iThink will import the data via the CSV file, and run the simulation. The results will appear in the graph in the model as shown below:

Updating the Data

In a real-world situation, the database that stores information we need to run a model is usually connected to another application that updates or changes the data. In the absence of a real-world application, I created a script that allows you to update the Simulation database and experiment with your own set of data.

In SQL Server Management Studio, open the script called UpdateTable.sql. Edit the values for each variable as desired.

USE [Simulation]
UPDATE Import_Data SET OnOrder=8, UnfilledOrders=4, TargetInventory=20, ActualOrders = 4
WHERE RowID=2
UPDATE Import_Data SET ActualOrders = 4
WHERE RowID=3
UPDATE Import_Data SET ActualOrders = 4
WHERE RowID=4

When you have finished editing the UpdateTable.sql script, execute the script to update the Import_Data table. Now you are ready to re-run the Transfer_Data.bat batch file and link the updated data with the Supply_Chain.itm model.

Summary

This example should give you a good idea of how the process of importing data from a database into an iThink or STELLA model can be set up and run automatically. Other database applications such as Oracle, MySQL and SQLite could also be used in a similar manner. I have not yet tried to set up examples, but these other databases also have command line interfaces for automating the transfer of data to CSV files that iThink and STELLA can connect to.

Stay tuned for part two of this post where I’ll provide an example of exporting model results to a SQL Server database.

The STELLA brand is targeted toward individuals in educational and research settings.  Supporting materials such as An Introduction to Systems Thinking with STELLA and sample models cover the natural and social sciences.

iThink, on the other hand, is targeted toward an audience of users in business settings.  An Introduction to Systems Thinking with iThink is written with the business user in mind and model examples apply the software to areas such as operations research, resource planning, and financial analysis.

Aside from the different program icons and other graphic design elements that go along with branding, there are just a few minor differences in the default settings for STELLA and iThink.  These differences are intended to pre-configure the software for the model author.  They do not limit you in any way from configuring the default setup to match your own individual preferences.

Below is a list of all the differences between the default settings for STELLA and iThink.

Opening Models

When opening a model with STELLA on Windows, by default, the software looks for files with a .STM extension.  Similarly, iThink looks for files with an .ITM extension.  If you want to open an iThink model using STELLA or vice-versa, you need to change the file type in the Open File dialog as shown below.

On Macs, the open dialog will show both iThink and STELLA models as valid files to open.

If you open a model with a file type associated with the different product than the one you are using, you’ll get a message similar to the one below warning you that the model will be opened as “Untitled”.  Simply click OK to continue.

Saving Models

When saving a model in STELLA, by default, the software saves the model with a .STM file extension.  Similarly, iThink saves model s with an .ITM extension.  If you’re using STELLA and want to save your model as an iThink file or vice-versa, use the Save As… menu option and select the appropriate type as shown below.

STELLA on Windows save dialog

STELLA on Mac save dialog

Run Specs

Since iThink is targeted toward business users who tend to measure performance monthly, the default Unit of time for iThink is set to Months.  It’s also easier to think about simulations starting in month 1 (rather than month zero) so we set the default simulation length in iThink to run from 1 to 13.  STELLA on the other hand, reports the Unit of time as “Time” and, by default, runs simulations from 0 to 12.

Run Spec Default Settings Comparison

Table Reporting

In a business context, financial results are generally reported at the end of a time period and the values are summed over the report interval.  For example, in a report showing 2010 revenues we would assume the values reflect total revenues at the end of the year.  In line with this assumption, the default Table settings in iThink include reporting Ending balances, Summed flow values, and a report interval of one time step.

In a research setting, scientists tend to prefer reporting precise values at a particular time.   For this reason, the default Table settings in STELLA are configured to report Beginning balances, Instantaneous flow values, and a report interval of Every DT.

Table Default Settings Comparison

STELLA or iThink

When choosing between STELLA or iThink, try to think about the kinds of models you intend to build and the problems you are looking to solve.  If your objective is to drive business improvement, chances are iThink will be a better fit.  If your purpose is to understand the dynamics of a natural environment or social system, STELLA will likely be your brand of choice.  Whatever you decide, both products will provide you with the exact same functionality and can easily be configured to suit your own preferences.

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.

Preparing the Model

The original model ran from 0 to 40 days. This needs to be changed to match the data, i.e., to 1 to 13 weeks. Changing this will force us to change the parameter values also, as described below.

The initial values of Susceptible and Infected must be changed to match the data set. The total number of people in the model should match the total population in the data. Summing the Hong Kong flu data set gives a population of 1035 people. The initial value for Infected is just the first reported value, 14. That leaves 1021 people that are Susceptible.

If you run the model using the original parameter values, you will find the infection is spreading far too quickly. To improve the chances of the calibration succeeding, it is necessary to adjust the parameters so that the results are in the ballpark of the historical data. To do this, it is useful to paste the historical data into a graphical function (called data in this model) and graph it alongside the model output. I first reduced the infection rate from 0.005 to 0.0015 and then adjusted recovery rate from 0.25 to 0.5 to get a closer fit to the data.

The model must also be modified to persistently import its parameters and export the time series we wish to fit to the historical data. The import sheet is most easily set up by creating an empty CSV file, exporting “all model variables” and “One set of values” to this file, and then removing the variables that are not needed from this file. The file can also be built manually. The final file, input.csv, looks like this:

infection rate,recovery rate
0.0015,0.5

You will also need to create an empty file called output.csv for the exported data.

Set up the import by choosing Edit->Import Data… and entering the settings shown below:

Add a table, called Table 1, to the model and add the stock Infected to the table. Set up the export by choosing Edit->Export Data… and entering the settings shown below. Note the “Use table settings” option is not used. This option cannot be used for calibration as it will restrict to precision of the output to match the table.

Save the model. It is ready to be calibrated.

Setting up PEST

PEST is a very powerful program. This post will only touch on the simplest capabilities that it provides. For more details and for more complex uses, consult the PEST manual.

PEST requires three text input files:

• The template file (.tpl), which tells it how to write the parameters to the model input file (input.csv in our case).
• The instruction file (.ins), which tells it how to read the data from the model output file (output.csv in our case).
• The control file (.pst), which tells PEST what to do.

Creating the Template File

This is the simplest file to create, as it only requires replacing data values in the model input file with markers so PEST can insert the parameter values it chooses.

Open the model input file, input.csv, in a text editor. Add the following line to the top of the file:

ptf #

This tells PEST this is a template file and that parameter fields, including the parameter name, will be surrounded by pound signs (#).

Replace each parameter value in the file with a short alphanumeric name for that parameter, surrounded by #s. Add any additional spaces that are necessary to make each field 13 characters long, including the pound signs. Using 13-character field widths, which is the maximum precision in single-precision floating point, is absolutely necessary to make certain the calibration succeeds (see end note about precision).

Save this file as input.tpl. The final file should look like this:

ptf #
infection rate,recovery rate
#infectrate #,#recoverrate#

Creating the Instruction File

The instruction file contains a set of instructions that PEST can follow to extract the model-generated data that corresponds to the historical data from the model output file. The model output file is a CSV file that has a header line followed by 13 lines, one for each week in the model run. Each of the 13 lines contains two numbers separated by a comma. The first value on each line is the week number and the second is the number of people who have been infected.

The contents of the instruction file, output.ins, appear below:

pif @
l2 @,@ !infected1!
l1 @,@ !infected2!
l1 @,@ !infected3!
l1 @,@ !infected4!
l1 @,@ !infected5!
l1 @,@ !infected6!
l1 @,@ !infected7!
l1 @,@ !infected8!
l1 @,@ !infected9!
l1 @,@ !infected10!
l1 @,@ !infected11!
l1 @,@ !infected12!
l1 @,@ !infected13!

The first line tells PEST that this is an instruction file and that the at-sign (@) will be used to delimit search keys. Each succeeding line tells PEST how to read the model output file, output.csv.

Each instruction line starts with a lowercase L, the PEST command to skip to the beginning of the Nth following line. The number one (1) is used to move to the start of the next line (interpreted as the start of the first line when at the beginning of the file), while larger numbers skip over an additional N – 1 lines. Thus, the “l2” in the first command, moves to the start of the second line in the file, skipping over the header. Each additional command uses “l1” to move to the start of the next line.

On every single line, it is necessary to skip to the second field (Infected), as we do not care about the week number. Surrounding the comma in at-signs (@,@) tells PEST to move to the next comma it finds (which will be the first comma on the line in all cases shown here). The exclamation points (!) delimit the name of the observation (historical data point) that corresponds to the data after the comma (called observation names in PEST). Note that each value of a time-series requires a unique observation name and that each observation name cannot be longer than 20 characters. In this case, there are 13 values, one for each week, in our historical data and our model-generated data, so the variables are named infected1 thru infected13.

In many cases, there will be far more data values that need to be assigned names. The instructions in this file, including the ever-incrementing names, can be quickly created in some text editors and also in Excel. The following formula can be pasted into cell A1 of Excel and filled down as far as necessary to create the proper instructions:

=CONCATENATE(“l1 @,@ !infected”, ROW(A1), “!”)

The above instruction names uses the name infected as the root of all observation names. Make sure to change this name in the above formula to match the variable name in your model. After you copy the text generated by Excel into your instruction file, make sure to change the first command so that it begins with “l2” rather than “l1”. If you fail to do this, PEST will not skip over the header in output.csv and will abort when it tries to read numeric values from the text that appears there.

There will also be cases where it is desirable to calibrate to multiple pieces of data. In this case, the commands in the instruction file will have to be modified to include the additional data. For example, if there is a second historical time series Treated and Treated appears after Infected in the model output file, each command would change to (using the fifth value as an example):

l1 @,@ !infected5! @,@ !treated5!

Creating the Control File

The control file contains a number of PEST parameters, the historical data, the command to run the model, and the names of the other files needed by PEST. It is a difficult file to create on your own. Luckily, PEST provides some command line tools to help. In particular, PESTGEN will generate the control file (with default settings) from two simpler text files:

• The parameter file (.par), which lists the input parameters to the model
• The observations file (.obf), which lists the historical time series; INSCHEK can generate this from the model output file and the instruction file Both files are fixed-format, which means the field widths and alignment must be properly respected. It is easiest to create the parameter file by editing an existing parameter file (PEST outputs a parameter file containing the final parameter values at the end of a calibration; the file used for this calibration is also included with the model).

The parameter file for this model, hkflu.par, appears below.

  single point
    infectrate    1.5000000E-03     1.000000         0.000000
   recoverrate    0.5000000         1.000000         0.000000    

The first line tells PEST to use single-precision floating point and to always include the decimal point when writing parameters to the model input file. The file then contains one line for each input parameter.

Each line starts with the name of the input parameter, as it appears in the template file, right-justified. This is followed by the initial value of that parameter, a scale factor (normally one) and an offset (normally zero). To add additional parameters, copy and paste the lines and make sure the columns line up with the existing lines.

The observations file contains all of the observation names with the historical data values for each. The command line tool INSCHEK, which comes with PEST, will generate this file by merging the observation names in the instruction file with data in a file formatted the same as the model output file. The instruction file has already been created. To complete this process, copy and paste the historical data into the second column (Infected) of the model output file, output.csv, and save it. Open the Windows command prompt window by choosing “Command Prompt” under All Programs->Accessories. Make sure to change the current directory to the directory containing your model and your PEST files. This is easily done by type “cd ” (the ending space is necessary) and then dragging the folder icon in the address field of the Windows Explorer window for that folder onto the command window. Reselect the command window and press Enter. The directory is now set. Now type:

inschek output.ins

Note that inschek will need to be preceded by the path to the PEST directory if you have not added the PEST directory to your PATH environment variable. This command will verify there are no syntax errors in your instruction file. If you get any errors, correct them and check again. Continue until there are no errors reported. Then type:

inschek output.ins output.csv

This will create the observations file, output.obf, from the data you pasted into output.csv.

It is now possible to create the control file, hkflu.pst. Type:

pestgen hkflu hkflu.par output.obf

If there are errors reported in the parameter file, fix them and try again. This command generates a default control file for case hkflu (short for Hong Kong flu – you can name it anything you want) with your parameter values and historical data. This file will require some editing before PEST can be run.

Editing the Control File

The generated file looks like this (to view, right-click on hkflu.pst, choose “Open With…”, and select WordPad):

pcf
* control data
restart  estimation
    2    13     2     0     1
    1     1 single point   1   0   0
  5.0   2.0   0.3  0.03    10
  3.0   3.0 0.001
  0.1
   30  0.01     3     3  0.01     3
    1     1     1
* parameter groups
infectrate   relative 0.01  0.0  switch  2.0 parabolic
recoverrate  relative 0.01  0.0  switch  2.0 parabolic
* parameter data
infectrate   none relative   1.500000E-03  -1.000000E+10   1.000000E+10
 infectrate     1.0000        0.0000      1
recoverrate  none relative   0.500000      -1.000000E+10   1.000000E+10
 recoverrate     1.0000        0.0000      1
* observation groups
obsgroup
* observation data
infected1      14.0000      1.0  obsgroup
infected2      28.0000      1.0  obsgroup
infected3      50.0000      1.0  obsgroup
infected4      66.0000      1.0  obsgroup
infected5      156.000      1.0  obsgroup
infected6      190.000      1.0  obsgroup
infected7      156.000      1.0  obsgroup
infected8      108.000      1.0  obsgroup
infected9      68.0000      1.0  obsgroup
infected10     77.0000      1.0  obsgroup
infected11     33.0000      1.0  obsgroup
infected12     65.0000      1.0  obsgroup
infected13     24.0000      1.0  obsgroup
* model command line
model
* model input/output
model.tpl  model.inp
model.ins  model.out
* prior information

This is also a fixed-format file, so you need to be careful when editing it. Luckily, only a few things need to be edited.

The “pcf” at the top of the file tells PEST this is a control file. The control data section sets a number of parameters that can stay the way they are unless the calibration does not work (in which case, you will need to dig deeply into the PEST manual to figure out what to change and how to change it). The parameter groups section, the observation groups section, and the observation data sections will not require any changes. However, the parameter data section needs to be adjusted and both the model command line and model input/output need to be specified.

For each parameter under parameter data, you must decide whether the changes made by PEST to the parameter will be relative (i.e., additive – the default), or factor (i.e., multiplicative). There are advantages and disadvantages to each, however, for a STELLA model, this simple guideline should help: change all time constants, rates, and multipliers to factor. In the SIR model, both variables are rates, so the word “relative” must be replaced with “factor ”. Note the two extra spaces after “factor”; the new text must completely fill the space formerly occupied by “relative”.

The default domain for each parameter is set to -1.0e10 to 1.0e10. You are likely to have better information (especially if it is a multiplier) and should enter it here in place of those numbers. I restricted the domain of infectrate to between 0.0001 and 0.01 and that of recoverrate to between 0.01 and 1.0. This was based on both my knowledge of these variables and some simple experiments I ran. The domains for these parameters could be much wider, as long as they stay above zero in both cases. This leads to the following changed lines (truncated after limits):

infectrate   none factor     1.500000E-03   1.000000E-04   1.000000E-02 …
recoverrate  none factor     0.250000       1.000000E-02   1.000000E+00 …

The model command line is changed to:

run.bat

The batch file run.bat contains commands to create the model output file and to run the STELLA SIR.stm model:

@echo off
echo > output.csv
start /min /wait “” “C:\Program Files\isee systems\STELLA 9.1.4\STELLA” -r SIR.stm

If STELLA was not installed to the default directory shown in the above command, run.bat will need to be edited to contain the correct path to STELLA.  The model output file needs to be recreated in the batch file (using echo) because PEST deletes the model output file before running the model.

Finally, the model input/output section was changed to name the template, model input, instruction, and model output files. Note these are not fixed-width fields.

input.tpl  input.csv
output.ins output.csv

Running PEST

Before running PEST, it is a good idea to verify that all of the input files are valid. This is accomplished with PESTCHEK:

pestchek hkflu

Fix any errors that are reported and repeat until all errors are gone. Note that PESTCHEK does not require the case name (hkflu) to be specified. Without the name, it will find all case names in the current directory (by searching for PEST control files) and check them all.

Finally, start the calibration process with:

pest hkflu

PEST will run until it detects no further changes or cannot proceed. Each iteration, it will report a value named phi. This is the sum of the squared residuals and should be decreasing each iteration. If PEST succeeds, the model input file will contain the best parameter values found and the model output file will contain the results obtained with those parameter values.

PEST outputs a number of files as well. As already discussed, the parameter file (.par) will be updated (or created) to contain the final successful parameter values. The run details file (.rec) will contain a log of the search, including parameter values tried and their relative success, as well as the final results compared to the historical data, with residuals. The 95% confidence interval for the parameters is also given. A number of other files are generated by PEST. See the PEST manual for further details.

The results of running PEST with this model and these initial parameters are shown in the graph below. The final parameter settings discovered by PEST were infection rate = 0.0015766646 and recovery rate = 0.82166505, with a final error (phi) of 4785.7 (compared to 773,059 with the original parameters).

Changing the Data Set

Shown below is the number of weekly deaths per thousand people in the UK due to the Spanish flu (H1N1) in the winter of 1918-1919.

To calibrate to this data set, it is necessary to change the initial values of the Susceptible and Infected stocks in the model. The historical data values in the PEST control file must also be changed.

The initial value of Infected is, again, the first data value (1 in this case). The initial value of Susceptible is, again, the remaining population (148). The model SIR2.stm contains these changed settings, as well as the new historical data set.

While it is relatively easy to regenerate and edit the PEST control file, it is just as easy to edit the few historical values. Control file spflu.pst contains the new dataset. It also uses run2.bat to run SIR2.stm instead of SIR.stm. Type:

pest spflu

Although it takes some time, the calibration almost succeeds using the same starting parameters as those for the Hong Kong flu. It fails in the end because the domain of infection rate is too narrow; infection rate gets stuck at 0.01. It is important to note that PEST reported success in this case. Looking at the graph in STELLA, however, made it clear that it only came close, while looking at the chosen parameters showed the boundary was reached. Changing the maximum value to 0.1 in spflu.pst fixed that problem. The final results are shown below.

This was generated with infection rate = 0.012101631 and recovery rate = 0.89054160. The error decreased from 2233.8 at the start to 110.6 with these parameters.

Troubleshooting PEST

Chapter five of the PEST User Manual includes a detailed troubleshooting section titled “If PEST Won’t Optimize.” If PEST is not changing parameters or not succeeding, this section will help you discover why.

Troubleshooting Data

In the examples used in this post, the model-generated data was one-to-one with the historical data. In many cases, data is not collected at a set frequency, does not coincide in time with model-generated data, and there are missing data points. While this was not covered in this example, PEST is capable of dealing with these cases. Missing data points can be handled by using !dum! as a placeholder in the instruction file (or additional lines can be skipped between points). PEST can also interpolate the model-generated data to match it in time to the historical data. Refer to the PEST manual for instructions to handle these cases.

End Note (about precision): While STELLA and iThink use double-precision for its calculations and PEST supports double-precision output, it is not possible to directly interface to PEST in double-precision mode as PEST outputs exponentials using “D” instead of “E” when in double-precision mode (as FORTRAN does). Such non-standard numeric formats cannot be read by STELLA. Double-precision can be used if necessary by adding a preprocessor to convert those “D”s in the model input file to “E”s before running STELLA and by adding a post-processor after running STELLA to change the “E”s output by STELLA to “D”s. Each parameter should also be given a field that is 23 characters wide, instead of 13, in the template file. Most applications will not need to use double-precision.

When a problem symptom appears, two options present themselves:  1) apply a short-term fix to the symptom, or 2) identify and apply a longer-term fix to the fundamental issue.  The second option is less attractive because it involves a greater time delay and probably additional cost before the problem symptom is relieved.  However, applying a short-term fix, as a result of relieving the problem symptoms sooner, reduces the desire to identify and apply a more permanent fix.  Often the short-term fix also induces a secondary unintended side-effect that further undermines any efforts to apply a long-term fix.  Note that the short-term fix only relieves the symptoms, it does not fix the problem.  Thus, the symptoms will eventually re-appear and have to be addressed again.

Classic examples of shifting the burden include:

• Making up lost time for homework by not sleeping (and then controlling lack of sleep with stimulants)
• Borrowing money to cover uncontrolled spending
• Feeling better through the use of drugs (dependency is the unintended side-effect)
• Taking pain relievers to address chronic pain rather than visiting your doctor to try to address the underlying problem
• Improving current sales by focusing on selling more product to existing customers rather than expanding the customer base
• Improving current sales by cannibalizing future sales through deep discounts
• Firefighting to solve business problems, e.g., slapping a low-quality – and untested – fix onto a product and shipping it out the door to placate a customer
• Repeatedly fixing new problems yourself rather than properly training your staff to fix the problems – this is a special form known as “shifting the burden to the intervener” where you are the intervener who is inadvertently eroding the capabilities and confidence of your staff (the unintended side-effect)
• Outsourcing core business competencies rather than building internal capacity (also shifting the burden to the intervener, in this case, to the outsource provider)
• Implementing government programs that increase the recipient’s dependency on the government, e.g., welfare programs that do not attempt to simultaneously address low unemployment or low wages (also shifting the burden to the intervener, in this case, to the government)

Moving from CLD to Model

The canonical CLD shown above can be confusing in two important ways:

1. The problem itself does not explicitly appear in the CLD, but is hidden in the left causal connection of balancing loop B2.  In fact, applying the fundamental solution reduces the problem, and therefore the symptom, rather than directly acting on the symptom.  In contrast, the symptomatic solution only affects the symptom, so the problem is never addressed and is free to generate more symptoms in the future.  This is shown in the revised CLD below.

   
   

   

2. The solution-induced side-effect affects the desire or ability to apply the fundamental solution, rather than the actual application of the fundamental solution as implied by the CLD.  This works in concert with the delay in balancing loop B2 to inform the decision about which solution to apply.

Even after acknowledging these two issues, the CLD remains difficult to translate into a model.  For example, what are the stocks and what are the flows?  I intentionally renamed what is normally called “symptomatic solution” and “fundamental solution” to “apply symptomatic solution” and “apply fundamental solution” to clarify these are flows, not stocks.  I likewise renamed “side-effect” to “solution-induced side-effect” to emphasize the side-effect was caused by applying the symptomatic solution.

Is the problem a stock or a flow?  It could be either, and that will determine whether the problem symptom is based on a stock or a flow.  In addition, since the problem is recurring, the fundamental solution will always correct a flow.  Finally, since the side-effect is continuously undermining the ability to apply the fundamental solution, i.e., it persists, it must be a stock.

This CLD hides one other detail that makes it hard to model:  the causal arrows leaving problem symptom imply a decision point; typically, only one path will be chosen.  The decision process is also implied in the CLD, being based on both the solution-induced side-effect and the delay, as well as a perhaps higher cost (not shown), to the fundamental solution.

It is unlikely the decision process will appear in a model constructed from this archetype.  When modeling the long-term impact of applying the symptomatic solution, there is no need to include either the fundamental solution or the decision process in the model.  In this case, the reinforcing loop R will guarantee the model continues to use the symptomatic solution until the system collapses (the Success to the Successful archetype again).  When modeling the long-term benefits to applying the fundamental solution, including short-term use of the symptomatic solution, the decision has also already been made and, again, the decision process remains outside the model.

In the end, I do not believe this is an archetype that is meant to be modeled in its entirety.  It allows discussion points to be raised and pieces of it can be seen in models.  Recognizing the pattern allows you to change the behavior and escape the potentially devastating results caused by the reinforcing loop.  In the following discussion, two models will be presented:  one that favors the symptomatic solution and another that favors the fundamental solution.  Only the second contains all of the loops in the CLD.  Neither one contains the decision process, as the decision was already made based on the purpose of the model.

Exploring Shifting the Burden with Credit Card Debt

People are often tempted to spend more than they earn, ignoring Mr. Macawber’s dire warning in Charles Dickens’ David Copperfield (with updated monetary units):  “Annual income $20,000, annual expenditure $19,500, result happiness.  Annual income $20,000, annual expenditure $20,500, result misery.”  Credit cards, unfortunately, make it very easy to achieve misery and most people do not see the danger of overextending themselves.

Consider someone who takes home $2000/month, but has managed to start spending $3000/month by a sudden change in lifestyle (new house, new car, etc.).  This person has $5000 remaining in savings at this point, no outstanding bills or debts, and a number of credit cards with varying limits over $5000.  It is clear that within five months, the savings will be used up.  What will the person do to cover their expenses?  One very common solution, the “symptomatic solution,” is to increase credit card balances by charging more purchases, including the use of credit card checks, and only paying the required minimum payment.

The CLD for this situation matches the Shifting the Burden archetype, as shown below with the symptomatic solution at the top and the fundamental solution at the bottom.

The model below (available by clicking here) implements this structure when the symptomatic solution is always chosen.

The model is divided into three sections, vertically.  The bottom section manages a bank account based on income and bill payments.  The center section manages spending, accruing bills, and the bill payment process.  The top section manages credit card debt (with no limits).

The model is annotated with all constant and initial values.  The balancing loop B1 corresponds to applying the symptomatic solution, i.e., drawing on credit to meet payment obligations.  The reinforcing loop R corresponds to the unintended side-effects created by this decision, i.e., mounting credit card debt that undermines the ability to apply the fundamental solution of balancing the monthly personal budget.  Note that the second balancing loop, B2, which corresponds to applying the fundamental solution does not appear in this model.  In fact, with regard to the solving the problem, this model is in open-loop:  there is no connection between spending and income (aka the right brain and the left).  While, ostensibly, the fundamental solution is to pay the bills from the bank account (there is a minor balancing loop there), there is no ability to do so without first bringing spending into line with income.  The real solution to the problem is to either reduce spending (shown in the CLD above), increase income, or do both.

A few assumptions allowed the model structure to be simplified without adversely impacting the outcome:

• The person has enough, or can obtain enough, credit to make it to the end of the simulation.  Thus, no credit limits, and no collapse, appear in this model.
• The system is continuous, i.e., there are no monthly payments per se, but payments are made throughout the month as one would expect from a normal stream of bills and a number of different credit cards all with different payment due dates.  The implication is that the person in question is faced with the problem symptoms (not enough money to pay their bills), and the choice between the symptomatic and the fundamental solution, on a daily basis.
• While there may be multiple credit cards, only one minimum payment cutoff is modeled.  This cutoff should be set to the sum of the minimum payments for all cards being used.

The model simulation produces the following results over three years.  Note how quickly the debt runs away, bringing the person to the brink of bankruptcy.

As expected, the bank balance is quickly depleted and the pipeline of bills increases to the monthly spending ($3000 plus the minimum credit card payment).  The credit card debt increases exponentially, and has reached over $35,000 after 30 months.  The monthly interest on this debt has risen to almost $600, nearly one-third of this person’s monthly income.

Investigating the Fundamental Solution

The fundamental solution (balancing loop B2) can be implemented by restricting the spending – including any credit card payment – based on the money available, as shown below (and available by clicking here).  In reality, it will take some time to adjust expenses downward, hence the converter time to react.  Note this is the inherent delay associated with applying the fundamental solution.  This means that it will still be necessary to borrow from the credit card, i.e., pursue the symptomatic solution in the short term.  This is a normal consequence of the Shifting the Burden structure.  The maximum credit card balance will then depend on how long it takes to adjust expenses down and how much extra money – beyond the minimum payment – is given to the credit card company.  There are, of course, lower interest alternatives to the credit card companies for short-term debt that would make it easier to reduce the debt balance.

The simulation results for the fundamental solution are shown below.  Note the credit card debt continues to grow even after spending has been controlled.  This is because the model continues to pay only the minimum payment ($45), which is less than the monthly interest charge (almost $90 by month 36).  To reduce the debt (almost $6000), an amount larger than the interest charge must be paid each month.

Forcing the minimum payment to be $100 (by increasing the minimum payment cutoff) is sufficient to both reduce the maximum debt to less than $5000 and to begin to lower the amount of debt.

Summary

Shifting the Burden is an important pattern that helps you recognize when you are caught in a never-ending cycle of applying a symptomatic solution.  As the first model shows, this can lead to long-term devastating results.  However, it is not true that symptomatic solutions should never be applied.  As the second model shows, sometimes they are needed as a short-term transitional fix to stabilize the system and/or buy time until the fundamental solution can be applied.  In these cases, a conscious decision has been made to sacrifice some long-term stability (increased debt burden over a long period of time in this case) in order to shore up the system until the long-term solution arrives.

The mid-term elections are now a month behind us and the political airwaves are still abuzz with commentary about the results.  Exit polls showed that unemployment was at the top of most voters’ list of issues, and that concerns about the federal government’s financial condition (record deficits and debt levels) were a hot topic as well.  Voters appeared to be asking “How can the federal government spend so much money and have so little positive impact on the nation’s economy?” 

The responses by politicians to such an important question are all over the map.  Democrats are claiming that economic conditions would have been much worse if not for massive federal bailouts and stimulus spending.  Republicans are touting the situation as a death knoll for the Obama platform in an effort to position themselves for 2012.  And the Tea Party movement has emerged to push for a roll-back of what they see as an intrusive and ineffective “Big Government”.

But, this political posturing reminds me that one of the true strengths of Systems Thinking is to force people to think very clearly and very operationally about the structure/behavior link embedded in such cases.  A little over a year ago, we sat down with Dr. Mark Paich, who used some very simple stock/flow language and some well-established principles of macroeconomics to lay out some relevant dynamics about the economic crisis and its aftermath:

• Why the collapse of the housing market made consumers re-evaluate their net asset position and hence started saving more of their incomes to pay off high interest credit card debt.
• How such actions on the part of consumers, in aggregate, kicked off a vicious cycle of decreased spending and contracting national output.
• Why government stimulus spending could close some, but not all, of the gap left by suddenly thrifty consumers, and that the recovery was likely to be a long, slow one.

We certainly don’t know how the future will play out, but the data suggest that consumers are indeed cutting back spending, and paying off debt.  (The Bureau of Economic Analysis has terrific historical data on household balance sheets and income.)  The unemployment numbers remain stubbornly high (around 9.5%), and although the recession is technically over, few economists are predicting rapid post-crisis economic expansion.

For a bit of clarity amidst all the rhetoric, you may want to check out Mark’s video offering.  His model and associated explanation do not provide a “magic bullet” of a solution, but they do provide some substance (and perhaps insight) to this vexing situation.  Now if only the politicians could follow suit!

To read a previous blog post about Modeling the Economic Crisis or view a 5-minute video trailer, click here.

In part one of this series I stated “A mental model is a model that is constructed and simulated within a conscious mind.”  A key part of this definition is that mental models are not static; they can be played forward or backward in your mind like a video player playing a movie.  But even better than a video player, a mental model can be simulated to various outcomes, many times over, by changing the assumptions.

Mental Simulation

Remember the example from part one of the child reaching for the hot stove?  One possible outcome we can simulate is that the child does not get burned.  We can simulate this outcome by altering our assumptions. We could include a parent in the room who rescues the child in the nick of time.  Or, we could simulate the child slipping just before reaching the stovetop because the hardwood floor appears slippery.  This kind of mental simulation allows us to evaluate what may happen, given different conditions, and inform our decision making.  We don’t have to make any decisions while looking at the picture, but imagine what actions you might take if the scene above was actually unfolding in front of you.

It seems effortless to mentally simulate these types of mental models.  Most of the time we are not even aware that we are doing it.  But other times, it becomes very obvious that our brain is working rather hard.  For example, looking at the chess board below, can you determine if the configuration is a checkmate?

It is indeed.  But I’ll bet it took noticeably more effort for you to mentally simulate the chess game than it did with the child near the stove scenarios.  Think about the mental effort that the players make trying to simulate the positions on the board just a few moves ahead in the game.

The paper “The Magical Number Seven, Plus or Minus Two: Some Limits on Our Capacity for Processing Information” by G.A. Miller (1956) established that people can generally hold seven objects (numbers, letters, words, etc.) simultaneously within their working memory.  Think of “working memory” as you would think of memory in a computer.  It’s like the amount of RAM we have available to perform computations within our mind.  And it’s not very much.  This means if people want to do any really complex information processing, they’ll need some help.  Over the last 50 years or so, the help has come from computers.  (In fact, IBM designed a computer specifically for playing chess, dubbed ‘Deep Blue’).

Digital computers have catapulted humankind’s ability to design, test, and build new technology to unbelievable levels in a relatively short period of time.  Space exploration, global telecommunication, and modern health care technology would not have been possible without the aid of computers.  We are able to perform the computation required to simulate complex systems using a computer instead of our minds.  Running simulations with a computer is faster and more reliable.

What makes a model useful?

Models that we can simulate using computers come in many forms.  For example, a model could be a financial model in a spreadsheet, an engineering design rendered with a CAD program, or a population dynamics model created with STELLA.  But what makes any of these models useful?  Is it the model’s results?  Its predictions?  I think the ability to explain the results is what makes a model truly useful.

Models are tools that can contribute to our understanding and decision making processes.  To make decisions, a person needs to have some understanding of the system the model represents.  A business finance model, for example, can be a useful tool if you understand how the business works.

Consider a model that does not provide any explanatory content, only results.  This type of model is often referred to as a black box.  It gives you all the answers, but you have no idea how it works.  People rarely trust these types of models and they are often not very useful for generating understanding.

The most useful models are structured so that the model itself will provide an explanatory framework that enables someone to ask useful questions of it.  Those questions may be answered by experimenting with the model (simulating) which, in turn, can help deepen a person’s understanding of the system.

This is an important feedback loop in a person’s learning process.  This feedback loop can be accelerated if the model provides explanations and can be simulated with a computer.

Transforming your mental models into visual models that are easier to understand and experiment with, will deepen your understanding, and help you communicate your models more effectively.

Modeling Dynamic Systems

Dynamic systems are notoriously difficult to understand.  A survey of opinions and policy proposals concerning the recent global economic crisis, health-care reform in the US, and climate change will yield a myriad of mental models about how these systems work, and what action we can take to improve them.  These systems are hard to understand because they challenge both our mental capacity to model, and to simulate them.  Systems problems are characterized by non-linearity, delays, and competing feedback loops, all of which are challenging phenomenon to understand.

As a result of the complexity that is inherent within dynamic systems, you’ll often find a lot of debate and mistrust of proposals to change them.  Think about how complex issues are presented in newspapers, blogs, on television and the workplace.  Typically, projections concerning what will happen, if we adopt one proposal or another, will be presented with little explanation of how the proposal will change the system’s behavior.  People will have difficulty engaging in meaningful discussion about how to actually solve problems without any common ground for understanding the dynamics of the system.

This is why developing a visual representation of how a system is structured underpins the philosophy behind Systems Thinking and the iThink and STELLA software.

Stocks and Flows

The mental model of the learning process pictured above is a visual representation of cause and effect relationships.  These cause and effect relationships need to be modeled operationally using stocks and flows to simulate them using STELLA.

A stock is an accumulation within a system.  Think of it as a bathtub.  Like a bathtub, stocks can be filled and drained.  In STELLA, stocks are indicated by a rectangle.  Stocks can be anything that accumulates—water, money, people, anger, and motivation are all examples of stocks.

A flow can add to, or subtract from, a stock.  Think of the faucet running water into the bathtub as it fills.  The stock is the accumulation of water in the bathtub and the flow is the rate at which the water is flowing in.

In the example below, the Population stock increases as people are born.   If people are born at a constant rate, then we’ll see a linear increase in the Population over time.

Similarly, the deaths flow pictured below, removes people from the population, causing the level to drop over time.

The clouds on the ends of the flows represent infinite sources or sinks.  In the model above, deaths flow into a cloud, which is outside the system boundary and therefore can accumulate infinitely.

Feedback Loops

Feedback loops are what enable systems to run themselves.  They are difficult to simulate mentally, but easy to construct using stocks and flows.  A connector (red arrow) indicates a dependency between two parts of a system.  In the model below, the births flow is dependent on the Population level and birth rate.

The birth rate is represented using a converter.  Converters can contain expressions or constants that can used to modify other parts of the model.  As you might suspect, the births flow is defined as Population * birth rate.  By placing the birth rate in a converter we can easily see that the flow of births is a function of both Population and birth rate.

This structure creates a reinforcing feedback loop.  Notice the population exhibits exponential growth when simulated over time.  There is nothing limiting the growth of the population.  One way to limit the growth is to add a balancing loop.

The relationship between the size of the Population and the death rate form a balancing loop.  The death rate changes based on the Population size and therefore limits the growth of the population.  If the population grows too dense, deaths occur at a higher rate.

Now that you have learned the basic building blocks (stock, flow, connector, converter) and understand how to map out feedback loops, we can build, and experiment with, a real model of a complex feedback system.

Predator-Prey Dynamics

A classic example of a complex feedback system is the dynamics exhibited by predator-prey populations.  A remarkable dataset from the Hudson Bay Company in Canada, of lynx and snowshoe hare pelt trading records, gives us a rare look at an isolated, natural system.  The records span almost a century, beginning in 1845.  From these records, we can infer the population levels of both species.

At first pass, it seems that there was little keeping the hare population in check other than the lynx population, (hares are the primary food source for lynx).  What is interesting here is the oscillation pattern.   How does the interaction between the two population systems generate this oscillation?  Let’s try to answer this question by constructing a model that we can easily understand and experiment with.

Understanding the System

What is your mental model of how a lynx population increases or declines over time?  My guess is that it may include a long-list of external factors (industrialization, weather, urban sprawl etc.) but the core of the mental model probably contains operational understanding; births and deaths ultimately control the population level.  Converting this mental model into an operational stock and flow model in STELLA looks like this:

The more lynx you have in a given population, the more lynx are born.  This is a reinforcing feedback loop.  As the population increases, lynx may start dying at a faster rate because they run out of some resource, such as food.  This is a balancing feedback loop.  Here are the results after I run the simulation:

We can see the population exhibits S-shaped growth.  The value at which the population levels off is the number of lynx that is sustainable within this particular system.

This simple model is congruent with most people’s mental models of population dynamics:  the population grows or declines based on the the limits of what is required to sustain it.  But let’s go a bit deeper and think about how this works.  What we have modeled is two feedback loops (reinforcing and balancing) competing with one another.

The graph above illustrates that the reinforcing feedback loop dominates at the beginning of the simulation, but as the population grows, the balancing loop takes control and dampens the growth.  Lynx are still being born when the population level stabilizes, it’s just that the number of deaths has become equal to the number of births.  Sometimes a balancing loop is referred to as a goal-seeking loop.  The level of the stock stabilizes when the goal is met.

Take a moment to think about how the feedback loops within the lynx and hare populations could interact with one another to generate the oscillation we saw.  We know that each population in isolation will grow until the carrying capacity of the system limits it.  What would happen if we modified the model diagram to show the dependency between the lynx and hare population?

I have connected the Lynx population with hare deaths and the stock of Hares with lynx death fraction via the hare density variable.  (The lynx death fraction changes based on the density of hares in the given area.)  I now have a useful model of a complex system: it is easy to understand how the model is structured, and I can experiment with the model by simulating different scenarios.

Experimenting

I have initialized the values in the model so that it begins in “dynamic equilibrium”—lynx and hares are dying and being born, but the population levels stay the same.  This is like a bathtub that has the faucet running at the same rate as the water leaving through the drain.  The level stays constant, but water is cycling through.  In this case, we have many generations of lynx and hares, but within the respective populations, the number of births and deaths are equal.

If we begin with the model in a steady-state, we’ll easily recognize change in behavior as we experiment.  Let’s explore the impact of trapping on the populations since we know both species were trapped for pelts when the data was collected.

I have added an additional outflow to the stock of Lynx that represents trapping.  I defined the outflow so that trapping begins after ten years.  At that point, 100 lynx are trapped and removed from the stock of Lynx.  Pictured below is the graph after the simulation is run with this new scenario:

And here is the graph after I set lynx trapped to 500 and run again:

In both experiments, the populations begin oscillating at the moment the trapping suddenly reduces the lynx population.  How can our knowledge of the feedback loops at play in this system help explain these results?

We know that in this model, the size of the hare population is limited only by the number of lynx.  When that number suddenly drops, the hare population increases because there are not as many lynx hunting the hares.  What happens to the lynx population?  Because there are more hares and less lynx, the lynx population increases as well.  (There is ample food available for the reduced lynx population.)  The reinforcing loop of births increases the population level, unconstrained by its balancing loop (Fig. 1).

Figure 1

Then something interesting happens:  the hare population peaks as the lynx population reaches its pre-trapping level.  At this point, the lynx population has increased to a level that limits the growth of the hare population (Fig. 2).  Lynx are eating hares faster than they can reproduce because of the increased hare density and number of lynx hunting.  In other words, the hare population’s balancing loop is limiting hare growth.

Figure 2

For a few years following the hare population’s peak, the lynx population continues to grow.  During this period, (as the hare population falls and lynx rises), the reinforcing births loop dominates the lynx population level (Fig 3).  There are ample hares available to prevent the lynx population from starving; the lynx balancing loop is not limiting their growth.

Figure 3

A few years later, the lynx population’s balancing loop kicks in as the number of hares plunges downward.  There are not enough hares left to support so many lynx, so they die at a much faster rate.  This rapid decline in lynx has the same effect on the hare population as trapping lynx: it begins another cycle of the populations overshooting the carrying capacity of the system, and subsequently collapsing.

How did our understanding of the model structure help us understand the results of our experiments?  We traced the cause of the oscillation back to the behavior of the feedback loops at play within the population systems.  In doing so, we gained a deeper understanding of the system that includes the behavior of two interdependent balancing loops.  This understanding can lead to more questions:  Will we observe the same behavior if we trap hares instead of lynx?  How can we stop the oscillation?  How will the oscillation change if we include the hare population’s food constraint in the model?

In the early 1900’s, the mathematicians Lotka and Volterra came up with a pair of non-linear, differential equations to describe the dynamics of a predator-prey system.  (They did not have the luxury of computers and software).  The same oscillating behavior is observed in economic systems and as a result the Lotka-Volterra equations played a role in the development of economic theory.  It is the understanding of the dynamics of oscillation that is useful here—much more so than the numbers.  This understanding enables us identify the structures that cause this behavior.  When we observe oscillatory behavior in the future, our mental models will be more accurate and we’ll be better equipped to make an informed analysis.

Better Mental Models, Simulated More Reliably

Our world is only going to become more complex and we need good thinkers helping to make it a better place for everyone.  Sharing systems models with people is a great way to increase their understanding of feedback loops and system dynamics.  Building models yourself is even more valuable, because you’ll learn a lot about your own mental models and improve them.  Practicing at this will help you to construct better mental models, simulate them more reliably, and communicate them more effectively.  Improving these skills will help inform what actions you take in the face of complex, systemic problems.

Connecting at the welcome reception

Earlier this month, we had the pleasure of hosting the 2010 isee User Conference in Providence, Rhode Island. During this amazing gathering of isee customers, partners, friends, and iThink/STELLA enthusiasts, we learned about the important work that is being done applying Systems Thinking to solve real-world problems, shared ideas, and made connections with one another.

For two and a half days, I saw participants immersed in keynote presentations, breakout sessions, and hands-on workshops.  As important, however, were the less structured round table discussions, poster presentations, and social activities.  Everywhere I looked, folks were engaged in conversation and connecting with one another.  Even better, I knew that many of those connections would continue to be fostered and developed well after the conference was over.

Presentations and models

Conference presentations and models are now available on the isee systems web site.  Please download and share these materials with your colleagues and friends. You’ll get a glimpse into the wide range of fields where Systems Thinking is being used to better understand the interconnections of dynamic systems including business, healthcare, education, energy, and the environment.  You can also download and listen to audio recordings of the keynote presentations:

• Business Tsunami: Surviving Disruptive Change in a Networked World
  by Henk Akkermans
• Education for an Interdependent World: Developing Systems Citizens
  by Peter Senge
• Learning and Teaching System Dynamics Modeling: A 20 Year Journey
  by Diana Fisher
• Generating Understanding of Social, Economic and Climate Dynamics in Interaction
  by Matthias Ruth
• The Youth Violence Systems Project
  by Steve Peterson and Paul Bothwell

A Conference Highlight

Steve Peterson describes the modeling process

One of the highlights of the conference for me was listening to the story that Steve Peterson and Paul Bothwell told about using Systems Thinking and dynamic modeling to help communities in high-violence Boston neighborhoods.

In their work with the Youth Violence Systems Project, there were two objectives:

• Improve understanding of community-based violence in Boston
• Help communities strategize and achieve sustained reductions in violence

What made this project different from other attempts to research and solve the youth violence problem in Boston was that it engaged “the community” in the development of the model.  From the get go, they included youth in the modeling process.  Gang members, in particular, turned out to be an important missing link to understanding violence and the dynamic system behavior.  Both Steve and Paul described some of the harsher realities of working with young people whose family members and friends were victims of violence.  The modeling process actually helped community members to articulate the “slippery slope dynamics” that move youth through the different stocks to gang involvement.  If you have a chance, I highly recommend listening to the audio recording!

Staying Connected

Participants engaged in round table discussions

Having time to interact with other participants was an important part of the conference experience.  It was wonderful to see the excitement and energy that is created when Systems Thinkers have an opportunity to connect with one another.  The cross fertilization that occurs so naturally between field experts, modelers, and educators was inspiring.

Please stay connected and let us know if there are other ways in which we can foster our growing community of STELLA and iThink modelers!

I first met Barry 20 years ago, when I had to trudge through the snow to get to his small wooden office.  I was building a discrete event model using STELLA and I wanted him to make some changes to the software so I didn’t have these “egg timer“ structures everywhere.  Barry was horrified with what I had done with his software and said words to the effect that it’s not meant for that, it was created to spread the gospel of System Dynamics.  Despite the fact that I was a civil engineer, he encouraged me to take a look at SD.  Such was his passion and conviction that 20 years later I’m still attending this conference.”

Tony Phuah accepts Scholarship Award from Mark Heffernan

Through most of his career Barry saw education as the key to spreading Systems Thinking.  As a teacher and a mentor, he dedicated much of his time to developing tools and methodologies for learning.  It is fitting therefore that this year’s award was presented to Tony Phuah, a Master’s student in System Dynamics at the University of Bergen.

Tony’s work includes an experimental study that explores the question: How can we improve people’s understanding of basic stock and flow behavior?  His experiment uses two different methods for teaching stock and flow behavior — the standard method (using graphical integration) and a method he calls “running total”.  Tony presented his paper at a parallel session during the conference and it can be downloaded by clicking here.  Although the results of his study favor traditional methods for teaching stock and flow behavior, we all should be encouraged by the work being done to try to improve Systems Thinking education and communication.  In Tony’s own words:

Speeding up ‘Systems Thinkers beget more Systems Thinkers’ growth will make us one step closer to Barry Richmond’s vision of a systems citizen world.”

Congratulations Tony and thank you Mark for helping us to celebrate Barry’s passion!

Applications for the 2011 Barry Richmond Scholarship Award will be available on the isee systems and System Dynamics Society web sites.  Check those sites for more information.

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.

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.

Running this model under Euler’s method leads to the behavior shown in the following graph.  As expected, Euler’s method tracks this function perfectly.  All three errors are essentially zero.

However, changing the integration method to Runge-Kutta 4 (known as RK4) does not give a good approximation (see graph below).  The relative error is 1.6%, the absolute error is $56 and the RMSE is $26, all rather large for RK4.  DT must be reduced a factor of 32, to 1/128, in order to bring the errors below 0.1% (0.001).

Continuous Compounding:  Runge-Kutta over Euler’s Method

Our friend Euler was a busy mathematician.  He is also famous for discovering the formula for continuously compounding interest.  If you increase the frequency of compounding from once a year to once a month to once a day (what most banks do now and called “daily compounding”) to once a second and so on until you are compounding every infinitesimally small instant (which should remind you of the dt from Calculus), the compound interest formula becomes:

Balance = Initial_Balance*e^(rate*time)

where e ≈ 2.718281828 is Euler’s number, rate is the annual fractional interest rate, and time is the number of years that have passed (a real number).  Using continuous compounding, our Balance at the end of year 3 is 1000*e^(0.1*3) ≈ $1350, $8 more than with annual compounding.

This function is known as continuous function, which, in layman’s terms, means its value smoothly changes with time, rather than experiencing jumps or gaps.  [Note the compound interest formula given in the first section is not continuous; since N must be an integer, this function only exists at integer values.]  While Euler’s method does not depend on a function being continuous, it does introduce an error proportional to DT when trying to integrate these equations (called truncation error).  Thus, the only way to reduce these errors is to reduce DT, requiring more computations.  Unfortunately, this only works up to a point.  Below that point, errors increase again due to round-off error in the computer’s finite representation of floating point values.

The graph below shows the effect of using Euler’s method to calculate continuously compounding interest.  Note it does not follow the actual curve, calculated with the EXP function, very well.  The relative error is 1.6%, the maximum error is $56, and the RMSE is $25.  If we cut DT in half, to 1/8, the errors also cut in half as we expect, to 0.8%, $29, and $13, respectively.  This will continue until DT reaches the limit of the machine’s representation, which is approximately 1/2048 for single-precision floating point numbers and 1/67,108,864 (1/2^26) for double-precision floating point numbers (note STELLA and iThink always use double-precision).  After this point, errors will increase again.  Note these values are theoretical; your mileage will vary.

There is a better solution than reducing DT.  The Runge-Kutta integration method assumes a continuous function and takes steps to follow such a function very closely.  Because it uses a weighted average of values across each DT-interval, it does not perform well with discontinuous functions, i.e., functions that have gaps or sharp changes, such as that produced with the STEP and PULSE functions (however, STELLA and iThink compensate for this so that these functions still perform correctly under Runge-Kutta).  Leaving DT at 0.25 and changing the integration method to RK4 produces the following graph, with the largest relative error being 0.0077%, the largest absolute error being $0.27 and the RMSE being $0.14.  To achieve this level of accuracy with Euler’s method, it is necessary to reduce DT to 1/1024.

The number after the RK is the order of the integration method.  Typically, but not always, higher-order methods will give smaller errors.  Euler’s method is a first-order method and RK4 is a fourth-order method.  Note, however, that Euler’s outperformed RK4 in the first example.  STELLA and iThink also provide a second-order method, Runge-Kutta 2 (RK2), which is a compromise between Euler’s and RK4; with the speed of today’s computers, it is no longer necessary to use and remains only for backward-compatibility.

Oscillations:  The Case for Runge-Kutta

Oscillations represent one area where RK4 really shines.  The following graph shows the results of using Euler’s method to integrate the sine function with DT = 0.25.  Notice it misses horribly.  The resulting function is clearly outside the range of the sine function [-1, 1].  The maximum absolute error is 0.25 and the RMSE is 0.15.  However, the maximum relative error is 340%!

Reducing DT by a factor of 128 to 1/512 helps somewhat (below).  The maximum absolute error reduces to 0.0020 and the RMSE to 0.0012, but the maximum relative error only decreases to 200%.

On the other hand, if we leave DT = 0.25 and change the integration method to RK4, it calculates the sine function quite well, with a maximum absolute error of only 0.0078 and an RMSE of 0.0045.  The maximum relative error is also acceptable at 2.2%.  Cutting DT in half to 1/8 reduces the maximum relative error an order of magnitude, to 0.34%.

Guidelines for DT

DT should initially be set to one-half of the smallest time constant in your model.  This minimal value, according to Nyquist’s Sampling Theorem, avoids undersampling your values.  [Undersampling means looking at your signal so infrequently that you distort the signal, or, more to the point, do not have enough information to reconstruct it.  For example, in the sine wave above, evaluating it every π will make you think you have a straight line at the origin, since all you will see are the zeroes.]

Finding the smallest time constant can be tricky.  They often appear in compounding or in draining processes.  For example, in the first two examples, the interest rate is the reciprocal of the time constant (it always appears in the divisor).  Thus, the time constant for this compounding formula is 10.  This suggests DT could start at 5, but there are other constraints.  In fact, we wish to compound every quarter, so DT must be set to 1/4.

Run the model with your chosen DT.  Then cut DT in half and try again.  If the results change, cut DT in half and try again.  If they do not change, the previously chosen DT is likely fine.

Whenever possible, choose a value of DT that is a power of 2 (1, 1/2, 1/4, 1/8, 1/16, … or 1, 0.5, 0.25, 0.125, 0.0625, …).  This is because most other fractions, for example 1/10, are repeating decimals in base-2 and so cannot be accurately represented on the computer.  Obviously if your problem requires you to use a DT of 1/7, in weekly simulations for example, do so, but understand that adding 1/7 seven times will not equal one.  There are also other cases where you might need to use a DT that is not a power of 2.  If your smallest time constant is 1/52, for example, choosing the closest power of 2 to (1/52)/2 = 1/104, which is 1/128, would be disastrous.  Your time constant would then be 1/64 sometimes and 1/32 other times, not at all close to what you need.  Make sure you can represent your smallest time constants as a multiple of whatever DT you choose.

Guidelines for Integration Method

Changing DT is often not sufficient.  Your system may exhibit instabilities or wild oscillations caused by the errors in Euler’s method.  It is therefore also wise to always run the model in Euler’s and then in RK4.  If there is a difference, stick with RK4 and proceed to test DT.

Euler’s method is noticeably faster than RK4 and works fine in many situations.  A system that has linear relationships and uses exponential functions that are not based on Euler’s number will work fine using Euler’s method.  For qualitative results, Euler’s will work well in a broad variety of situations.  Euler’s is also the preferred method when you are using discontinuous or discrete functions (but note the discrete functions in STELLA and iThink are Runge-Kutta-aware and should work reasonably well under RK4).

Any continuous system that is based on exponential functions using Euler’s number as a base, and this includes trigonometric functions, oscillations, and logistic functions, will be better served by RK4.  In these cases, you should generally use this method from the start.

There is one additional type of system you may run into that neither of these methods will solve very well:  stiff systems.  In simple terms, a stiff system is one that has time constants of wildly different magnitudes, for example, 0.005 and 10,000.  The smaller time constant suggests a DT of 0.0025 which is far too small for the part of the system with a time constant of 10,000.  It is best to reframe your problem to avoid these situations when possible.  For example, when simulating such a system, the part of the system with a time constant of 10,000 can often be considered constant over a shorter period of interest.  Over long periods of interest, it is often possible to aggregate the behavior of the part of the system with the shorter time constant, using averages instead of exact values.

Note:  The integration method and DT are both set in the Run Specs, accessible from the Run menu.  Additional information can be found in the chapter on DT in the Technical Documentation.

It’s been a little over 10 weeks since the Deepwater Horizon oil rig explosion that has resulted in a constant flow of oil into the Gulf of Mexico.  Oil is now beginning to impact the economy of the Florida coast.  Some estimate that the disaster could cost nearly 200,000 tourism jobs.   Efforts  to remove the impacts on the environment, including a massive rescue of manatees could cost billions.  The ability to truly restore the environment to pre-April 20 conditions is beyond that of mere mortals. It is an event truly unprecedented.

Who’s to Blame?

The current media focus has centered on an activity I refer to as the Find the Knucklehead Game. The idea is that there must be someone out there to blame, tar and feather, and that if we just remove the idiot from the system we’ll never have this problem again. Finé. Complete. Case closed! (And there was great rejoicing…)

However, if we begin to apply the systems thinking paradigm, there’s another analysis that might suggest we are in a long-term trend that is still just ramping up – and that if we don’t take action soon the impact on the economy and environment may be much worse.

A Systems Thinking Perspective

Two of the skills required to practice systems thinking are Dynamic Thinking and 10,000 Meter Thinking. If we look at the history of oil extraction from the balcony, as a long-term trend over time, we see something that may be quite useful. Let’s apply a concept developed by Charles Hall (SUNY-College of Environmental Science and Forestry) called EROI (“energy return on investment”).  EROI is the ratio between the energy we receive – to run our transportation system, heat our buildings, run electricity generators – and the amount of energy required to get the raw material out of the ground and process it into usable form. In the early 20th Century, oil was easy to extract. In many places it was just below the surface in the wide open fields. The EROI was 100:1 in 1930 – 100 units of energy received for 1 unit of energy extracting/processing. Since then there has been a marked decline, and as the United States passed peak oil production (the maximum production rate) in the early 1970s, and as we’ve begun importing most of our oil, the EROI for oil in the US is approximately 20:1.

Charles Hall's Energy Return on Investment Graph

Why the marked decline? Because the purest, easiest to extract/process oil has likely been found and burned (leading to too much in our atmosphere, but that’s another point). We now need to go to more distant places (no longer fields) – offshore, into the Arctic – to find oil. Often the oil is no longer in purest form; converting the “oil goo” from the tar sands in Alberta, Canada is a prime example.

Continuation of a Trend

Using the systems thinking perspective, we may conclude that the “event” of Deepwater Horizon is an inevitable continuation of a trend (that includes the Exxon Valdez). And that the trend is likely to get worse as our need (thirst) for oil increases – which will occur as developing nations continue trying to catch up to the developed world, and as we hope for “economic recovery.”

This leads to several questions in my mind.

1. If it’s true that the easiest to extract and process oil has been found and used, what’s this indicate about the risks we’ll incur as it requires increasing effort to get oil in the future? More deepwater drilling? Environmental degradation/damage to get oil like tar sands extracted?
2. How much oil really is left – how long will we continue to have this cheap resource?
3. Further, if we are facing increasing risks and costs for every unit of oil we acquire, when will that begin to have severe impacts on our economy (as the per unit cost of the stuff we love – iPods, cars, flying to Europe – increases, the profit margins of our globalized corporations will decrease)?
4. And when will we decide that we need to develop a way of living that is independent of this resource? To recover from our “addiction to oil.”

How Much Oil Is Left?

I’ve developed and published an isee NetSim model that you may use to explore questions regarding how long the resource will last and the implications of the economy on that length. You may explore it here.

Ultimately, we need to stop the Find the Knucklehead Game and instead recognize that it is we – in the collective sense – that are responsible. Not just for the recent disaster, but the long-term trend and its consequences (including the unintended consequences of climate change). The system is behaving as it is designed. It is up to us to design a different system!

Barry was fully engaged in bringing systems thinking to everyone.  He saw how this powerful way of thinking could help people to better understand society’s most pressing issues and make the world a better place.  Barry saw K-12 education as one of the keys to creating a better world.  He spoke often about educating young people to become ”systems citizens” and preparing students for the complex problems they would have to face.  Much of his time was devoted to training teachers to incorporate systems thinking into curricula and pedagogy.

A couple of years ago, Barry’s daughter, Joy Richmond, began spearheading an effort to create a book in honor of her father.  Joy invited a group of us together to talk about some ideas for the book and come up with a plan to make it happen. The first idea we discussed was writing the book that Barry himself had intended to write.  Barry left plenty of notes and even had a working title for a book about systems thinking called Traces.   We all agreed that it would be much too daunting to try to write a book for Barry, so we decided to have a book written in tribute to Barry by friends and colleagues who share his passion for systems thinking.

Steve Peterson, Corey Peck and Khalid Saeed were all part of that original discussion and eager to contribute by writing a chapter.  Each had a story to tell about using Systems Thinking in their work and why it is so important in an increasingly interdependent world.  What better way to honor Barry than writing a book that helped get the word out about systems thinking!

Shaping the Book

Lees Stuntz, Executive Director of the Creative Learning Exchange, was also in on the discussion and excited about asking educators influenced by Barry to contribute their stories. Before we invited other authors however, we wanted to provide some guidelines that would tie the book together and give it a more meaningful context.  I think it was Steve who came up with the idea to use the critical thinking skills first outlined in an article Barry wrote for the System Dynamics Review titled “Systems Thinking: Critical Thinking Skills for the 1990s and Beyond”.   We agreed the systems thinking skills would provide a good foundation for the book and each author could then choose a few of the thinking skills to emphasize when telling their story.

Countless hours of writing, editing, and designing later, Tracing Connections: Voices of Systems Thinkers was born.  Published in partnership with the Creative Learning Exchange, proceeds from the book will fund scholarships that offer learning opportunities for educators to use systems thinking and system dynamics in K-12 education.  The response so far has been excellent and we are pleased to be funding scholarships to help educators attend the ST/DM Conference later this month.

A Chapter for Everyone

What is especially nice about the book, is that you don’t need to read each chapter in sequence.  Since the authors’ experiences range from education and research to business and public policy, there’s sure to be a chapter for everyone.  Click on the link below to view the chapter by Frank Draper titled “Teaching by Wondering Around: Learning About the World Naturally”.  Frank tells a wonderful story about how Systems Thinking has transformed the way he teaches science to high school students.  After reading it, you’re going to wish you could enroll in one of Frank’s field science classes at Catalina Foothills school district in Tuscon, Arizona.

Teaching by Wondering Around by Frank Draper

Animal Temperature Model

Table of Contents with full list of chapter titles and authors

For more information or to order a copy of Tracing Connections, visit http://www.iseesystems.com/tc

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)

Forio Simulate—Like YouTube for Simulations

Since the debut of the isee NetSim software we have collaborated with Forio Business Simulations to provide free hosting for models published online.  Last year, I saw a demo of Forio’s new simulation service, Forio Simulate, in the early development stage.  It was awesome.  I couldn’t believe it—they added all of the Web 2.0 features I’d been dying for, and a slick, fun to use interface.  I knew we had to make isee NetSim work with it.

The new service is social by nature—it features keyword tagging, ratings, comments and built-in model sharing.  It’s like YouTube for simulations.  You can browse simulations and games that others have created, find similar ones by keyword tags, and embed simulations on blogs or web pages.  So you don’t have to be actively publishing your own simulations, or even be a modeler to get a lot of value from the site.  In fact, it’s a great way to introduce people to modeling and simulation.  Just post a link!

And if you are a modeler, you can get even more out of the service.  Most of the simulations published to Forio Simulate allow you to download or copy the underlying model that powers it.  Here at isee, we’ve uploaded the original iThink and STELLA models used to create all of our sample isee NetSim simulations.  That way anyone can run the simulation online in a browser and download the source model to dig into if they want more.  Take a look at isee systems’ models published on Simulate by clicking on the screen shot below.

The simulation overview page allows anyone to download the source model, rate, comment, add tags and embed the simulation on other web pages.

Anyone can use free hosting on Forio Simulate if the model contains 400 or fewer equations and is shared with other users of the service.  If you want to protect the model, resell a simulation, use a custom URL, etc., you can sign up for a premium hosting plan.  It’s easy to select a hosting plan when you upload your model with isee NetSim.

Graphical Input Device

Use the GID to change the price assumption

One of the most requested features for isee NetSim has been the graphical input device, or “GID” as we like to call it.  GIDs offer a simple way to allow people running your simulation to change the model assumptions.  Often it’s easier to sketch a curve of the behavior of a variable than to input numbers directly.

For example, the model below is a generic renewable-resource model that enables you to experiment with yield and price of a renewable-resource via GIDs.  The model is from our “Thinking in Systems” online course.  I published the model to Forio Simulate and then embedded it here.  Click the “Review model” link to take a brief tour of the model.

You can change the price elasticity and “yield per unit capital” curve by first clicking on one of the GIDs.  A window will pop open that displays the details of the curve.  Click and drag within the grid to draw a new curve.  Press “Ok” to apply the changes.  Take a few moments to play around with the model below.

[Editor's note:  If you cannot see the embedded model below because you are reading this post in a feed reader or email, please visit the post page here]

What I find interesting after playing with this model is what happens when you start to push the yield up.  At some point the efficiency improves to the point where the renewable-resource gets completely wiped out.  Sometimes the inefficiencies of harvesting a resource can be a good thing—it gives the resource time to “catch up”.  This may be a useful point to consider when thinking about sustainable systems.  What happens when you make the resource inexpensive?

A GID for the web

You may have noticed that the GID within isee NetSim does not allow a user to type numbers to define data points or change the x and y ranges.  We decided this simpler version of the GID was a good fit for models published to the web and would cover most use cases.  The GID allows users to easily change the assumptions by sketching a behavior pattern and not really concern themselves with the data.

If you’ve been dying to publish models with GIDs online, now is your chance.  While you’re at it, check out Forio Simulate and start contributing to the growing simulation community.  You can download a trial of isee NetSim 1.1 here.  If you already own isee NetSim, visit your “My Software” page on the isee systems website to download the upgrade.

Background

The April 26, 2010 article in the New York Times titled “We have met the enemy and he is PowerPoint” has created quite a stir. It is particularly telling that three days after its publication, it is the most emailed article on their website! The most interesting aspect of the article to me, as a system dynamics practitioner, is the publication of a system dynamics map on the US Counter-insurgency strategy as the example (i.e. visual sound bite) demonstrating why PowerPoint is so problematic. This is actually a poor example of the author’s point, since it is not PowerPoint, and because the map was shown out of context.

Although the diagram doesn’t portray how or why PowerPoint is misused, it does demonstrate some reasons why system dynamics maps and models are not more broadly used to communicate systemic issues. In this post, I will describe what issues a PowerPoint paradigm creates and how system dynamics can address those issues; more importantly, I will show why the STELLA and iThink software have features such as storytelling and web publishing in order to help people develop deeper, more systemic understanding of the complex problems humanity must address.

The Problem with PowerPoint

I don’t have anything against using PowerPoint. Those of you who have taken one of my webinars for isee systems know that I rely heavily on the software in my instruction and facilitation. I think there are inherent software limitations that combine with a cultural paradigm, that lead to its misuse. Currently, I see it promotes the following approaches to problem solving:

1. Narrow focus in space and time – due to limited screen real estate
2. Passive absorption of information of data – lazy learning, not experiential
3. Simplistic bullet point thinking – linear thinking, focusing on factors in a non-operational way

This all creates confusion between reducing complication and simplifying complexity. The world is a dynamically complex place, and thank goodness for that! Picture the blandness of a world that is simple, where everyone thought and acted the same, where you always knew exactly what would happen because it was so simple. Boring! On the other hand, dynamic complexity makes it difficult to resolve what currently appear to be intractable problems, such as environmental degradation, poverty, global economic turmoil. Living in a dynamically complex world necessitates finding ways to simplify complexity to its essence, making manageable and useful mental models.

That’s why people are drawn to lists (e.g., bullet point and linear thinking), believing it simplifies complexity; just give me a list of what’s wrong or what to do! What lists do well is remove complication, but they also remove the dynamic essence of reality, often making mental models that are less than useful.

System dynamics

System dynamics is an approach to building understanding that expands boundaries, looks at the world as comprised of feedback loops, uses a visual language that promotes operational thinking, and creates active learning. It’s a terrific approach to counter the many problems inherent in applying PowerPoint paradigm!

All of the above helps develop useful mental models that are both simplified and still capture the essence of reality. However, taking a map out of context – even one much simpler than shown in the article – and including it in PowerPoint will not create understanding, only confusion! When I’m in front of a group and have enough time, I will always draw it up on a flipchart or board, to bring the group along with its unfolding. The rapid feedback creates an engaged group capable of learning. But in the absence of time, or if you need to communicate to people “on their own time” you will find features in STELLA and iThink invaluable!

I’ve published a map to the web with the isee NetSim software to demonstrate how you can use system dynamics to create online experiential learning labs. Take a tour of the map below to see how the stock/flow language and Storytelling can overcome the passive absorption of bulletized information that PowerPoint facilitates.

Click on the image below to make sense of the map

Another interesting perspective from Linda Booth Sweeney on the New York Times article can be found on her blog.

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.

Eliminating Cross-Sector Connections

To convert this sector-based model to modules, it is first necessary to remove these cross-sector connections.  This is done with the ghost tool.  The stock Population needs to be ghosted and placed in the Resources Sector and the stock Natural Resources needs to be ghosted and placed in the Population Sector.  The end result is shown below.

Creating the Modules

We are now ready to create, in a new model, one module for each sector, using the Module tool.  The result is shown below.  Although you could manually draw the connections between these modules, we will let the software create the connections for us.

Each module needs to be associated with its own model file.  This is done by double-clicking on each module and choosing “Create Model…”.  Double-click Population, choose “Create Model…”, and name the model file “Population”.  Repeat for Natural Resources, but name the model file “Natural Resources”.

Pasting Structure into the Modules

It is now a simple matter to paste the structure from each sector into each module.  Using the selection tool in the sector frame, select all of the structure in the Population Sector and copy it by choosing Edit->Copy.  Open the Population module in the new model and paste in the structure by choosing Edit->Paste.  The module contents will appear as shown below.  Note that Natural Resources is no longer a ghost; there is no connection to the other module.

Repeat this process with the Resources Sector and the Natural Resources module.  The module contents will appear as shown below.  Note that Population is no longer a ghost.

Connecting the Modules

All that remains is to connect the modules.  Inside the Population module, the stock Natural Resources needs to be replaced with a ghost of the real stock.  Inside the Natural Resources stock, the Population stock must be replaced with a ghost.

There are two ways to do this:  manually or semi-automatically with the Ghost tool.  I tend to use the Ghost tool most of the time.  [The documentation chapter “Working with Module Inputs and Outputs” (available by clicking here) explains how to forge these connections manually.]

To connect the modules, first go into the Natural Resources module and select the Natural Resources stock with the Ghost tool.  Go up one level and then down into the Population module.  The cursor should still be a stock.  Click next to the existing Natural Resources stock to deposit the ghost.  Run a connector from the ghost to resources\person (as the stock does) and then delete the stock.  Drag the ghost over to approximately where the stock was.  Open resources\person and press “OK”.  An error will appear and “Natural_Resources” will be selected.  This happens because Natural Resources is no longer connected to resources\person.  Click on the ghost’s name (Natural Resources.Natural Resources) in the “Required Inputs” list and press “OK”.  The Population module is now properly connected to the Natural Resources module.

Repeat this process to connect a ghost of the real Population stock to the Natural Resources module.  After finishing, the Population module contents appears as shown below.  Note that Population has a double border to signify it is providing its value to another module.

The contents of the Natural Resources module appears as shown below.  Note the double border around the Natural Resources stock, signifying that it is providing its value to another module.

Finally, returning to the top-level, we see that the software automatically created the required connections between the modules.  Below, the connections have been embellished with their polarity to make the balancing loop obvious.  [The final model is available by clicking here.]

A Cautionary Tale

This mechanical process for converting a sector-based model to a module-based model is very appealing due to its simplicity.  It works very well with models that have a small number of sectors (up to about five), but will fall apart beyond that.

Last summer, I helped Khalid Saeed (Worcestor Polytechnic Institute) convert a sector-based model of Jay Forrester’s Urban Dynamics model to modules.  A direct translation using the above process led to the following top-level diagram:

It should be clear that everything is connected to almost everything else.  There is no way to rearrange this diagram to clear up this mess.  For a short time, I flirted with the idea of module ghosts as a solution, but it does not take too long to see that ghosts cannot solve this problem either (indeed, they make it worse).

Recently, I decided to follow my own advice on module hierarchy (see this post).  I also remembered the famous article “The Magical Number Seven, Plus or Minus Two: Some Limits on Our Capacity for Processing Information” by G.A. Miller (1956).  It explains that we cannot keep more than five to nine things simultaneously in our heads, depending on the person.  So I may be able to absorb the meaning of a diagram with seven modules, but someone else may only be able to understand one with five modules.  For this reason, we should strive to keep the number of modules in a diagram below seven (which also may help explain why the above method works well for a small number of sectors).

Rearranging the Urban Dynamics model into three levels, I was able to simplify the top-level of the diagram to the following, which shows all of the important high-level feedback loops:

There are also lower-level feedback loops within each of these modules.  The one within workforce is just as complicated:

The model logic rests within these second-level modules.  For example, the module underemployed mobility contains the following logic:

In writing and teaching people about Systems Thinking, we often refer to “mental models”.  For some people, this comes as a bit of a surprise, because the context usually involves building models with the iThink or STELLA software.  They don’t expect us to start talking metaphysically about thinking.  “Is this about philosophy or modeling software?” they may wonder.  The software is actually a tool to help construct, simulate and communicate mental models.

Let’s define the term model: A model is an abstraction or simplification of a system.  Models can assume many different forms – from a model volcano in a high school science fair to a sophisticated astrophysical model simulated using a supercomputer.  Models are simplified representations of a part of reality that we want to learn more about.  George Box stated: “Essentially, all models are wrong, but some are useful”.  They are wrong because they are simplifications and they can be useful because we can learn from them.

So, what is a “mental model”?  A mental model is a model that is constructed and simulated within a conscious mind.  To be “conscious” is to be aware of the world around you and yourself in relation to the world.  Let’s take a moment to think about how this process works operationally.

Imagine that you are standing outside, looking at a tree.  What happens?  The lenses in your eyes focus light photons onto the retinas.  The photosensitive cells in your retinas respond by sending neural impulses to your brain.  Your brain processes these signals and forms an image of the tree inside your mind.

So at this point, we’ve only addressed the mechanisms by which you perceive the tree.  We have not addressed understanding what a tree is or considered changes over time.   We are dealing with visual information only.  There is nothing within this information that tells you what a tree actually is.

What makes the image of a tree in your minds click as an actual tree that exists right there in front of you?  This is where mental models kick in and you start to think about the tree.  The tree is actually a concept of something that exists in physical reality.  The “tree concept” is a model.  Understanding the concept of a tree requires more information than is available through sensory experience alone.  It’s built on past experiences and knowledge.

A tree is a plant.  It is a living thing that grows and changes appearance over time, often with the seasons.  Trees have root systems.  Trees use leaves for photosynthesis.  Wood comes from trees.  I can state these facts confidently because I have memories and knowledge of trees contained within my mental models.  Mental models contain knowledge and help us create new knowledge.

Take a look at these images for a few moments and then think about what is happening inside your mind as you look at them.

My guess is that with each image you thought about what happens next.  If so, you were actually simulating a mental model forward in time.  The images don’t show you what happened next, but you can probably make a pretty good guess.  This guess is the result of simulating a mental model of what is depicted.  And, you can simulate different outcomes.  In the image on the right, I like to simulate a parent rushing into the scene and scooping up the child before he gets burned.

When we think about the boy catching the ball, the blocks knocking one another over, worrying about the toddler getting burned, we are applying knowledge that we hold in our minds to help simulate a mental model of what is depicted in the image.  This is all going on subconsciously, so we’re really not aware of it as it is happening.

Thinking About Systems

The human mind is very good at simulating mental models of our immediate physical reality.  Things get harder when we start thinking about abstract systems.

A market is a good example of an abstract system.  In a market system, price acts as a signal of aggregate demand for a commodity.  You can’t “see” a market like you can “see” a tree in front of you.  A market does not exists in a particular physical location.  A market is an abstract concept that exists in the collective minds of all who participate in it.  Even though a market does not exist physically, markets have an enormous impact on our lives nonetheless.

When the global economic crisis hit in late 2008, retailers began to struggle financially because consumer purchases declined rapidly.  People were worried about the economy and started saving money instead of spending it.  This started happening just before the holiday shopping season – a make-or-break period for many retailers.  So, in an effort to boost demand, retailers began dropping prices.

This led to price deflation because consumers saw prices dropping rapidly and began delaying purchases as a result.  The outcome of simulating their mental model of the market informed their decision making: “I should wait to buy this because the price keeps dropping”.

This mental model paints a pretty picture for consumers over the short-term:  low prices in a down economy.  As the deflationary dynamics play out over the long term, the picture becomes bleak.  As prices spin downward, profits decline and businesses are forced to layoff workers or close up shop entirely.  As unemployment increases, consumers perception of the stability of the economy decreases and they spend even less.

Economists and policy makers use sophisticated computer models to help them understand markets.  On the other hand, consumers use simple mental models when making purchasing decisions.  The more sophisticated models inform policy makers of the long term consequences of consumers cutting spending, so they react by trying to jump-start spending with stimulus programs.  In the US we’ve seen a few of these programs during 2009: the “Cash for Clunkers” rebate program, the First-Time Home Buyer Tax Credit, and the Social Security payroll tax cut.

Often, it is hard for us to define the optimal boundaries for a mental model.  We tend to have a narrow focus and act on short-term dynamics within our mental models.  For example, in the model above, our understanding changes when we expand the boundaries to include profits and layoffs.

However, we are not generally very good at mentally simulating complex systems with interdependencies, lots of variables, and delays.  This is where software steps in.  Using Systems Thinking software, we can transform our mental models into operational models that we can simulate more reliably using a computer.  This not only helps us create new knowledge and understanding, but helps us construct better mental models in the future.

We’ll explore the role of software in constructing, simulating and communicating mental models in part two of this series.

For more on using modeling and simulation to understand the economic crisis, see the “Understanding the Economic Crisis” video presentation.

With this in mind, we collaborated with our consulting and training partner, Lexidyne LLC, to create a new series of video-based presentations focused on common but often misunderstood problems that can be conceptualized, expanded, and then explored using Systems Thinking.  We recently released the first video in this series — Understanding the Economic Crisis presented by Dr. Mark Paich.

Judging from its title, you might think Understanding the Economic Crisis presents a huge complex model of the macro economy.  To the contrary, dynamic modeling expert, Mark Paich, begins with a very simple model of something we all can relate to — the individual consumer.

Stock and flow map of an individual consumer's balance sheet

Mark expands upon the model and shows how a sudden drop in housing prices affects individual consumption.   As you might expect, when Total Net Worth falls, the individual responds by spending less.  When housing prices fall, home equity loans no longer provide the purchasing power for big ticket items like cars, vacation homes and big screen TVs.

The real surprise however, comes when Mark further expands the individual consumer model to include the economy as a whole.  When everyone’s net worth decreases at the same time a phenomena known as the “Paradox of Thrift” occurs. The paradox states that if everyone tries to save money during times of recession, total savings for the overall economy may fall.  The dynamics generated by adding elements of the macro economy to the model are indeed surprising.

Mark’s easily understood model leads to some real insights concerning the policy implications for an economic recovery.   It also provides a great example of how Systems Thinking can be used to deepen your understanding of a complex issue in order to make better decisions.  If you haven’t seen the video, I highly recommend it.  The following trailer highlights some of the key points in Mark’s presentation and will give you a taste of the full presentation.

(If you cannot see the video below in your RSS reader, please visit the post page)

For more information or to purchase Understanding the Economic crisis, click here.

Last time, we explored the effects of Type 1 rivalry (rivalry between different companies in a developing market) on the Bass diffusion model by replicating the model structure.  This part will generalize this structure and add Type 2 rivalry (customers switching between brands).

Bass Diffusion with Type 1 Rivalry

To model the general case of an emerging market with multiple competitors, we can return to the original single company case and use arrays to add additional companies.  In this case, everything except Potential Customers needs to be arrayed, as shown below (and available by clicking here).

For this example, three companies will be competing for the pool of Potential Customers.  Each array has one-dimension, named Company, and that dimension has three elements, named A, B, and C, one for each company.  Although each different parameter, wom multiplier, fraction gained per $K, and marketing spend in $K, can be separately specified for each company, all three companies use the same values initially.  All three companies, however, do not enter the market at the same time.  Company A enters the market at the start of the simulation, company B enters six months later, and company C enters six months after that.

Recall that the marketing spend is the trigger for a company to start gaining customers.  Thus, the staggered market entrance can be modeled with the following equation for marketing spend in $K:

STEP(10, STARTTIME + (ARRAYIDX() – 1)*6)

The STEP function is used to start the marketing spend for each company at the desired time.  The ARRAYIDX function returns the integer index of the array element, so it will be 1 for company A, 2 for company B, and 3 for company C.  Thus, the offsets from the start of the simulation for the launch of each company’s marketing campaign are 0, 6, and 12, respectively.

This leads to the following behavior:

Note that under these circumstances, the first company to enter the market retains a leadership position.  However, companies B and C could anticipate this and market more strongly.  What if company B spent 50% more and company C spent 100% more than company A on marketing that is similarly effective?  This could be modeling by once again changing the equation for marketing spend in $K, this time to:

STEP(10 + (ARRAYIDX() – 1)*5, STARTTIME + (ARRAYIDX() – 1)*6)

The effect in this case is quite different.  All three companies end up with approximately equal market shares.

Bass Diffusion with Type 1 and Type 2 Rivalry

Once a market is established, companies will also make an effort to win competitor’s customers away.  A previous post showed how to model switching between brands using arrays.  Adding that structure to this model gives a model (available by clicking here) that implements both type 1 and type 2 rivalry.

The specific mechanism that convinces customers to switch brands is not explicitly modeled here so that the structure matches that in the earlier post.  Instead, this model uses constant switching rates; the earliest entrant, company A is given the most favorable values.  However, to make the model more dynamic and also useful, there should be factors that influence the relative attractiveness of the two brands, for example, marketing spend, product quality, product features, and customer service.

Giving an advantage to the early entrant to the market has the expected effect of boosting their customer base at the expense of the others.  This another example of Success to the Successful.

Note the scale of this graph had to be changed relative to the earlier graphs to show all of company A’s customers.  The effect of the customer switching is so strong that even increasing marketing for company’s B and C does not change the end result.  Of course, the switching factors in this model are constant and increased marketing should indeed change them.

This does, however, show that it is important to intervene in the correct place to reach your desired outcome.  Once the market is mature, the dynamics are dominated by the switching effects, not the process of gaining customers.  Thus, efforts to change market share by changing the rate of gaining customers are ineffective.  Energy should instead be focused on reducing the number of customers who switch to other brands and increasing the number of customers gained from other brands.

This model does not model the decline of the market, through customer loss, since we are mostly concerned with growth and maturation.  If desired, the loss would be modeled with an outflow from Customers as shown in the earlier brand switching post.

Finally, this model assumes that Customers belong to one company alone, i.e., no one uses the product from more than one company at one time.  Sharing customers is known as type 3 rivalry and is beyond the scope of this post.

A Note About Revenues

Gaining customers is, of course, not the entire story.  The purpose of gaining customers is to create revenue for your business.  There are two opportunities to generate revenue in this model:

1. From your existing customer base.  This stock-driven revenue stream is typical for manufacturers of non-durable goods such as bleach and razors, or for subscription-based or service businesses, such as magazines, cable TV, or beauty salons.  If your customer base is large, your revenues will be high.  The Bass diffusion model shows a growing customer base and therefore a growing revenue stream (see curve 1 below).
2. When a person is converted from a Potential Customer to a Customer or a customer switches from one brand to another.  This flow-driven revenue stream is the only source of revenue for manufacturers of durable goods, such as cars, computers, and refrigerators.  However, it often also applies to software companies.  Since the flow rate determines the revenue stream, if the rate at which you gain customers is large, your revenue will be high.  Note, however, that the Bass diffusion model shows us this rate falls sharply as the market matures, as the number of Customers saturates (see curve 2 below).  To keep revenues high in a mature market for this kind of business, it is necessary to entice customers away from competitors, build in product obsolescence (so your customers have to buy again), or enter new markets.

A Final Aside on the Bass Diffusion Model

In his excellent book Strategic Modelling and Business Dynamics, John Morecroft explains (on page 173) that the effect of marketing in the Bass diffusion model has a stronger purpose than just to start the ball rolling.  Without the marketing effect, market growth through word-of-mouth alone is much slower than any business would desire, or typically experience.  The following graph compares the growth with and without the marketing effect (Customers is initialized to one in both cases).

Without the marketing expenditure, the conversion of Potential Customers to Customers does not peak until the month 91, whereas, with marketing expenditure, it peaks in month 50, in almost half the time.  This is also reflected in the time it takes the market to mature:

I think it is safe to say that we all were very impressed by the quality of the students’ work and how well they understood the dynamics associated with the issues they were presenting.  Perhaps more striking however, was seeing how empowering modeling real-world issues is for young people and the enthusiasm they share for their work.

Now everyone can see the effect that modeling real-world issues has on students at the CC Modeling Systems web site. Dedicated to helping educators bring dynamic modeling into the classroom, the web site features videos of students presenting their work as compelling evidence to the value of incorporating System Thinking and system dynamics into curriculum.

You’ll be amazed to see what 14-18 year olds are capable of and the excitement they exuberate when addressing challenges such as:

• How much do carbon emissions need to be reduced in order to stop global warming?
• Can automakers cope with the increasing demand for hybrid gasoline-electric vehicles?
• What is the impact of introducing an invasive species into an ecosystem?
• Will China have sufficient water resources to sustain population and agricultural growth?
• How can we improve the blood donation process in the United States?

Students are eager to understand the world better and are more than capable of building and understanding relatively sophisticated models in their attempts to understand the dynamics of real-world systems.

—Diana Fisher

Educators and administrators considering dynamic modeling curricula typically face challenges. No matter how compelling the evidence that Systems Thinking and the system dynamics methodology engages students and takes them to a higher level of reasoning, it is still difficult to justify without tying it to National Standards.

The CC Modeling Systems web site devotes an entire section to detailing very specific 21st Century Skills and National Standards addressed by curriculum that incorporates building system dynamics models.   Much of the homework has been done aligning this type of work to standards in the following subject areas:

• 21st Century Skills
• Economics
• Health
• Math
• Science
• Social Studies
• Sustainability
• Technology

Many thanks to Diana Fisher for sharing her students and her experiences teaching dynamic modeling with all of us!

To learn more about the modeling course that Diana teaches, I recommend the following links:
Recorded webinar presentation by Diana Fisher
Modeling Dyamic Systems: Lessons for a First Course

In part one of this series, I explained the Limits to Growth archetype and gave examples in epidemiology and ecology. This part introduces the Bass diffusion model, an effective way to implement the capture of customers in a developing market. This is also used to implement what Kim Warren calls Type 1 rivalry in his book Strategy Management Dynamics, that is, rivalry between multiple companies in an emerging market.

The Bass Diffusion Model

The Bass diffusion model is very similar to the SIR model shown in part one. Since we do not usually track customers who have “recovered” from using our product, the model only has two stocks, corresponding loosely to the Susceptible and Infected stocks. New customers are acquired through contact with existing customers, just as an infection spreads, but in this context this is called word of mouth (wom). This is, however, not sufficient to spread the news of a good product, so the Bass diffusion model also includes a constant rate of customer acquisition through advertising. This is shown below (and can be downloaded by clicking here).

The feedback loops B1 and R are the same as the balancing and reinforcing loops between Susceptible and Infected in the SIR model. Instead of an infection rate, there is a wom multiplier which is the product of the Bass diffusion model’s contact rate and the adoption rate. If you are examining policies related to these variables, it would be important to separate them out in the model.

The additional feedback loop, B2, starts the ball rolling and helps a steady stream of customers come in the door. If you examine the SIR model closely, you will see that the initial value of Infected is one. If no one is infected, the disease cannot spread. Likewise, if no one is a customer, there is no one to tell others how great the product is so they want to become customers also. By advertising, awareness of the product is created in the market and some people will become customers without having encountered other customers who are happy with the product.

The behavior of this model is shown below. Note it is not different in character from the SIR model or the simple population model.

Implicit Type 1 Rivalry

In the face of competition for an emerging market, multiple companies will be fighting for these same potential customers. Ostensibly, they all follow the Bass diffusion model. However, if the exact dynamics of the competition is not important to the purpose of the model, it is possible to simply model the reverse pressure the competition is placing on your company’s growth.

In this version of the model (which can be downloaded by clicking here), both the wom multiplier and the marketing effectiveness (fraction gained per $K) are restricted by competitive pressure. The variable competitive pressure is an S-shaped graphical function meant to mimic the behavior of competitors in a developing market. In the beginning, there is no competition. As the market grows, more and more competitors enter the market until it becomes saturated. It is a simple time-series function, as shown below, but could instead be tied to the success of your company.

The variable competition is just 1 – competitive_pressure, i.e., it is a multiplier that reduces the effectiveness of wom and marketing spending as competitive pressure grows. Finally, a switch called competition switch has been added. When this is set to zero (the default), the model behaves as it did previously, i.e., without competition. When it is set to one, competition is active and the behavior changes to that shown below.

Explicit Type 1 Rivalry

Finally, another competitor can be explicitly modeled by replicating the structure (or by using arrays). This model is shown below (and can be downloaded by clicking here).

In this structure, each competitor can be given their own explicit parameters, which can also vary as the simulation progresses based on market share or product attractiveness. When the rival spends half as much on marketing as your company does, your company wins more market share:

In this simulation, the rival enters the market at the same time as your company. This is usually not the case in the real world. If all parameters are the same for both companies, but the rival company enters the market one year later, the end result is identical to that shown above. The transient response, however, is necessarily different:

This model does not include the process of developing a market for a product nor does it model persons leaving the market. I.e., it is assumed that the total market size is predefined and never changes. For a completely new product, the potential market will likely grow, so we may want to include this effect. However, when modeling a developing market, it is unlikely we are interested in the shrinking of the market that occurs after maturation.

The model also does not model losing customers to competitors (Type 2 rivalry) or customers who use products from multiple competitors (Type 3 rivalry). An example of modeling Type 2 rivalry is shown in a previous post.

My next post in this series will expand the above model to include Type 2 rivalry and to use arrays.

Rolf Olsen, Workshop Participant

I was very excited about a last-minute chance to attend the introductory iThink/STELLA workshop, but to be honest, on the flight to Colorado Springs, I started to become apprehensive.  Who was I trying to kid?  Sure, I’d heard the terms “stock” and “flow” and I understood their roles as the nouns and verbs of the software.  I’d even read a few chapters in Barry Richmond’s Introduction to Systems Thinking.  But the first time I started up the software and stared at that blank workspace, I had no clue where to begin!  Adding to my anguish, I was quite certain there would be others there who were much smarter than me and really knew what they were doing.

In college I spent most of my time and energy studying English and French, language, literature, cinema, art history, and so forth. I managed to avoid all higher math like the plague (although I did reasonably well in basic statistics).  My engineer father often reminded me that my degree in Humanities prepared me for almost nothing.  After college, I stumbled into a career in marketing – quite fertile territory for exploring system dynamics and modeling, as it turns out.  I spent a few formative years in an ad agency and at a regional banking system, before finding my stride marketing and managing nonprofit arts and culture organizations. Today I work in marketing and communication in a large academic medical center.

For years I’ve used spreadsheets to model various ‘what if’ scenarios.  In the arts, I used spreadsheets to create budgets and set ticket prices, always seeking ways to better predict revenue from ticket sales at different prices, for different types of performances (e.g., modern dance, string quartet, jazz ensemble), or on different days of the week.

Preparing for the iThink/STELLA workshop, I decided I’d like to try to model demand in a market area for laser vision correction surgery, popularly known as LASIK or PRK.  That seemed simple enough.  I might be able to bluff my way through this workshop after all!

In the end, all my concerns and fears of inadequacy were for naught.  Corey Peck is a remarkable teacher.  Within the first hour of the first day, he’d already explained that almost everybody who starts iThink for the first time takes a deep gulp at the sight of that blank workspace, wondering how to begin.  That was comforting, to say the least.

Others in this workshop included a nurse from Nigeria, two men from a large consulting firm that counts the military among its principal clients, three people with separate interests in systems used to deliver electrical energy, two from the IS planning department of a large national insurance company, a woman from an Oregon health care system, and two women from Hong Kong who teach English at a university there.  And there was me, the health care marketing geek.

Despite the diversity of our backgrounds and fields of endeavor, somehow, over two and a half days of workshop, Corey always managed to guide and instruct us, as a group and individually, in a way that was neither too complex and esoteric nor too simple.  He obviously comprehends dynamic systems so well, and knows the software so thoroughly, that no matter what type of model one of us presented, Corey could quickly offer clear and helpful guidance.  I never felt lost, and I always found myself learning something as I eavesdropped on Corey’s conversations with others in the class.

Our other instructor, Mark Paich, was a protégé and colleague of some of the pioneers of Systems Thinking, including Barry Richmond and Peter Senge.  Mark’s tales from years of consulting work with clients, such as General Motors, were fascinating and, for me, they helped place the concepts we were exploring within a much larger context.  Although the scope and complexity of some of the projects Mark described were daunting, these stories helped me appreciate more clearly what is possible when you begin to view almost any situation you encounter as a dynamic system.  And the stories were fascinating and fun, providing an excellent “brain break” from the rigors of our lessons with Corey!

By the end of the first day, we’d all constructed our first working models, and begun to sketch out the models we hoped to build.  As I mentioned, I wanted to try to create a simple model to explore dynamic demand, based on several variable factors, for LASIK vision correction surgery.  I’d start with the population of our market area and then determine the number of people who wear glasses or contacts.

• Of those people, how many are aware that our hospital offers LASIK surgery?
• How many of the “aware” group might be willing to attend a free information session and preliminary evaluation?
• Of all qualified candidates, how many would be willing and able to actually go through with the surgery?
• Then, after successful LASIK surgery, how many people would tell their friends, neighbors or colleagues how happy they were, thus spreading awareness and influencing a new group of candidates to have LASIK surgery?
• And what might be the impact on surgery volume of marketing efforts at various levels of expenditure?  How about a discount or referral incentive?

I started out pretty simply because, as Corey reminded us often, “There’s plenty of time to make it more complex later.”  That proved to be great advice.  And in the end, I went home with a rudimentary model that seemed to work.  Well, mostly.  It needs more work and I need more practice.  And when I mentioned to Corey that I hoped to figure out how to factor in the “word of mouth” effect into my model, he promised to send me information about the “Bass Diffusion Model” that addresses that very topic.  And sure enough, an email arrived a few days after I’d returned home.

In all, my experience at the introductory iThink/STELLA workshop was excellent, even transformative.  I now often find myself thinking in terms of stocks, flows, converters and connectors. I don’t literally model every situation I encounter, of course, but I see applications all around me. I know I have an enormous amount yet to learn, the introductory course covered just the basics for me. But now I have more confidence and I feel inspired to move forward, because I’ve had a tantalizing glimpse of what’s possible with a proper grasp of key concepts, the right software tool – and remarkable teachers.

As the first anniversary of the isee Blog approaches, we thought it would be interesting for folks to see the list of our most popular blog posts.

Top Ten Posts of 2009

1. Modeling H1N1 Flu Outbreak
2. Modeling Customers Switching Between Brands
3. Modeling a Watershed with Arrays
4. Matrix Arithmetic
5. Spatial Modeling with isee Spatial Map
6. “Thinking in Systems” book inspires online course
7. Physics Textbook 2.0
8. Insight-based Model Investigates the Housing Crisis
9. Building a Health Care Model Hierarchically
10. C02 in the Atmosphere Behaves Like a Bathtub

The Limits to Growth Systems Archetype, also known as Limits to Success, combines growth with an exogenous or endogenous limit.  This Systems Archetype was formally identified in Appendix 2 of The Fifth Discipline by Peter Senge (1990), but made its first prominent appearance in World Dynamics by Jay Forrester (1971) and then The Limits to Growth by Meadows, Meadows, Randers, and Behrens (1972).  The Causal Loop Diagram (CLD) is shown below.

Real growth processes have inherent limits to growth.  Identifying these limits can help avoid problems in the future, whether the problem is overpopulation, increasing demand for a product that cannot be met, or growing a business in a mature market.  When growth is desired, but limited, it is always better to find ways to increase the limit before pushing for more growth.  Excessive growth in the face of a limit often leads to collapse.  Driving the system to the point of collapse can erode the ability to continue after the collapse, for example, by reducing the production capability of a piece of farmland or destroying the reputation of a company.

Classic examples of limits to growth include:

• The collapse of the deer population on the Kaibab plateau and on St. Matthew Island due to overpopulation and the attendant overgrazing of their habitat
• The overshoot and collapse of the human population on Easter Island
• Overgrazing in the Sahel region of Africa by cattle herders
• Overfishing of the oceans by fishermen
• The collapse of People Express due to sharp customer growth combined with slow personnel growth
• The sharp exodus of America Online subscribers after an intense marketing campaign increased the number of subscribers far beyond their capacity
• The contraction of the world economy in 2008 due to limiting oil supplies
• The productivity of staff deteriorating as a company grows, due to increased interactions and reporting overhead
• Business growth limited by the size of the potential market
• Yeast cells in the fermentation process, who suffer from both the loss of exogenously supplied sugar and the increase of endogenously produced pollution

Exploring Limits to Growth with the SIR Model

The SIR (Susceptible-Infected-Recovered) model, also known as the Kermack-McKendrick model, is used to predict the spread of an infectious disease.  The disease is ultimately limited by the number of potential hosts.

In this model (available by clicking here), the Susceptible population becomes infected based on their rate of contact with an Infected population and the probability of being infected (infection rate combines these parameters).  After some amount of time (1/recovery rate), the disease runs its course and people either recover or die, moving into the Recovered population.

The infection feeds on susceptible people.  As the number of susceptible people falls, the number who can become infected also falls.  Since infected people eventually recover, the number of infected people will ultimately fall.  This behavior is reflected in the graph below.

The sharp increase in Infected, followed by a dramatic drop is known as overshoot and collapse.  It is a characteristic behavior of the Limits to Growth Archetype.  The stock of Recovered people is also limited – by the total population.  It exhibits a different characteristic behavior of Limits to Growth:  S-shaped growth, so called because the graph of this growth pattern looks like the letter S.  This is also known as logistic growth.  In the beginning, growth is very rapid, but as the limit is approached, growth slows and then stops.  In feedback terms (refer back to the CLD), the rapid growth on the left side of the S occurs because the reinforcing loop dominates while the right side slowdown occurs because the balancing loop dominates.

There are two separate limits to growth:  one on the sustainable level of the condition of interest and one on its rate of growth.  The former is defined by a system’s carrying capacity (named limiting condition in the CLD).  As Khalid Saeed, Professor of Social Science and Policy Studies at Worcestor Polytechnic Institute, recently reminded me, the carrying capacity of a system is always a flow.  Any system that requires something, either energy or mass, to persist cannot exist beyond the level at which that something can be provided.  The second limit, that on the rate of growth, is defined by the system slack, which is always a stock.

In the SIR model, the infectious disease feeds off of the susceptible people.  The stock of Susceptible people is therefore the slack in the system.  The growth of the disease is limited by the contents of this stock.  However, the carrying capacity of the system is zero, as there is no inflow to Susceptible to provide an endless supply of victims for the disease.  The net result is that eventually, and after a sharp rate of growth, the number of infected people will also drop to zero.

In either a closed or a quarantined community, this may be true.  However, ease of travel creates an inflow to Susceptible, i.e., a carrying capacity, and thus the disease can persist.  The graph below shows what happens when an inflow is added Susceptible and set to 10 people per day (10% of the initial population).

Rather than going to zero, both Susceptible and Infected wind up at non-zero steady-state values that reflect the new carrying capacity (10).  Note in particular that this value is not 10 for either stock; it is about 50.  Note also that the Infected peak has risen from roughly 17 to 90 due to the inflow of susceptible people.

For more information about the spread of infectious disease, read our blog post on Modeling H1N1 Flu Outbreak.

Exploring Limits to Growth in Population Dynamics

Another example occurs in population dynamics.  The simple example shown below demonstrates the impact of a limit (click here to download this model).

The food available is a simple linear relationship that reduces available food as population increases.  The death rate increases as the available food per member of the population falls.  The behavior of this model is shown below.

While this shows the desired behavior, an exogenous limit (food available) enforced as the population grows, this model has no slack and no explicit carrying capacity.  The lack of an explicit stock to model the slack in the system will result in an incorrect transient response (i.e., the system’s response to a change in its variables) when something happens that might require that slack, for example, a drought or a surge in the population due to immigration.

A more accurate and useful, though more complex, model appears below (and is available by clicking here).

Both the carrying capacity (replenishment) and the slack (Food Available) are explicitly modeled.  The behavior is quite different:

There is now a proper transient response to an increase in the birth rate.  The slack, Food Available, fuels growth in the Population until the slack is exhausted (relative to the carrying capacity), at which point the population collapses to its approximate equilibrium value – the value supported by the carrying capacity.

Next time, I will explore limits to growth in the context of business development models.

We’ve been working with Mark over the past year to incorporate dynamic modeling and computer-based simulation into the STEM curriculum at Ocean Lakes High School in Virginia Beach.  Ocean Lakes is serving as a demonstration school for a broader initiative in Virginia to use Modeling and Simulation (MODSIM) as an engine for 21st Century workforce development.

One of our goals for the pre-conference workshop was to provide participants with STELLA models and lessons that they could use immediately in their classrooms.  To attain this we put together a packet of sample lessons that would cover the spectrum of STEM education courses—Science, Technology, Engineering and Mathematics.

Science: Exploring an H1N1 flu outbreak model

Educators at all levels have found that physical activities or games can often be a good way to introduce students to STELLA models. Activities can be a lot of fun and provide a physical model for kids to make the connection with the more abstract computer-based simulation.

In our workshop, we decided to introduce a STELLA model of the H1N1 flu outbreak by first engaging the class in a simple exercise that demonstrates the spread of infectious disease using cups of water.  Each student was given a cup of water and a pipette.  One cup of water however, was contaminated with sodium hydroxide.   We then began the process of walking around and dropping a pipette filled with liquid from our own cup into the cups of fellow students.  At the end of several rounds, Mark (our teacher) put a drop of phenolphthalein solution into everyone’s cup.  If the liquid turned pink, we were infected.

With a new understanding of how quickly infection can spread, the class was ready to be introduced to a STELLA model of the H1N1 flu outbreak.  The core structure of the model is based on the same model that epidemiologists use to track a population through the various stages of infection including Susceptible, Exposed, Infected and Recovered (SEIR).

The simulation provided several controls for students to experiment with different policy options for controlling the spread of the virus.   Along with the model, we provided a sample lesson with a set of questions to help guide students through their exploration.  To download the H1N1 model and sample lesson, click here.

Technology: Inquiry-based learning with STELLA simulations

It was fascinating to watch how each student in our workshop experimented differently with the model and began asking their own set of questions.  STELLA simulations and computer technology provide the perfect platform for students to learn using an inquiry-based approach.  Rather than being told how something works, students can discover for themselves by exploring “what-if” scenarios and finding out what happens. For teachers, this can mean shifting from a traditional teacher role to that of a facilitator or coach, guiding students in their inquiry without knowing in advance the path they’ll choose.

Engineering/Physics: The Pendulum Story

Another STELLA model that we introduced to workshop participants was one that Mark and I modified for the Physics Flexbook project earlier this year.  The original version was developed by Diana Fisher in her book Lessons in Mathematics: A Dynamic Approach. The model provides a practice field for learning how the variables of simple harmonic motion are related.  Controls in the simulation allow you to explore what effect, if any, string length, initial displacement, and pendulum ball mass have on the amplitude, period, and frequency of the pendulum’s motion.  To download the pendulum model and sample lesson, click here.

Mathematics: Algebra Word Problems

You may remember from your own experience in math class that word problems typically give students a lot of trouble!  It is especially difficult to understand what is important in the problem and how to translate words into mathematical equations.  STELLA can help by allowing students to create a visual representation of the problem and making it easier to understand the symbolic representation.

For the workshop, we used the following word problem as an example of a typical Algebra I assignment:

Imagine your class is going to try to raise $400 by making school T-shirts. There is a $150 set-up charge for the T-shirt design that you have designed. Once the design is set, it costs $4 for each T-shirt. You feel it is possible to charge $10 for each T-shirt. How many T-shirts do you have to sell before you break even, i.e., make enough money to cover your costs?

Together we created a STELLA model that provided a visual representation of the word problem and gave everyone some hands-on experience building models from scratch.  It also inspired some participants to think about ways they could expand or change the model to answer additional questions that weren’t part of their immediate assignment. 

To download the T-shirt model and lesson, click here.

We’d love to expand our library of STELLA lessons and models so that we can share them with teachers.  If you’ve got any ideas of things you’d like to try in your classroom, please feel free to contact us – we’re here to help!

The STELLA model is based on the SEIR compartmental model that epidemiologists use to model the progress of an epidemic.  SEIR models divide the population into compartments: Susceptible, Exposed, Infected and Recovered.  These ‘compartments’ translate nicely into stocks within the STELLA model where we can observe the dynamics of the spreading virus.

While developing the model we decided to explore some strategies that schools are pursuing to limit the virus’ spread.  We wanted to know if the “stay at home” (when you are sick) policy would be effective in the case where vaccines are not available quickly enough, (which as of November 2009 is the case).

Take a look:

Click the ‘Simulate’ link on the home screen above and try some different scenarios.  Be sure to click the ‘How does this simulation work?’ link for a guided tour of the model behind the simulation.

As you experiment with the simulation, consider the following:

• How does varying “% vaccinated” effect the number of sick students?
• How many days do infected students need to stay home to have a significant impact on the spread of the virus within the school?
• What impact does the “% effectiveness of vaccine” have on the flu outbreak?
• What combination of decisions results in the lowest number of sick students?  Are these decisions realistic in a real-world setting?

Note: Each time you dial in parameters and press run, a new plot will be added to the graph so you can compare the effectiveness of the different decisions.  Clicking on the blue reset button will clear the graph and reset all Knobs to their default value.

If you think this simple model is useful, feel free to share it or embed it on your own website; just click the sharing icon in the lower right corner.  If you want to dig deeper into the STELLA model you can download the model by clicking here.  You can open the model with STELLA 9.1, or the free isee Player.

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.

Given these probabilities, we can find the number of customers switching out with the equation:

switching_out_probability[Products]*Customers[Products]

We can then use the exact same formulation as previously for switching to:

switching_out[Products]*switching_in_probability[Products, Products]

and switching in:

ARRAYSUM(switching_to[*, Products])

Note that all customers that switch between two brands also pass through a stock called Customers Switching.  This was Steve Peterson’s idea:  If the equations in the model are correct and all rows of switching in probability add up to one, as they ought to, this stock will always be zero.  If it is ever non-zero, you immediately know there is something wrong with the numbers in switching in probability.  [Note:  Due to truncation errors during the simulation, the stock will not identically equal zero, so it has been globally scaled from –1 to +1 using the Range Specs.]

Expanding the Model

The structure given above models rivalry for existing customers quite well.  However, it does not include flows for new customers entering the market or existing customers leaving the market.  The full structure for acquiring and retaining customers looks like this (and is available by clicking here):

New customers are gained by each brand based on their respective market shares.  Each brand also has a fraction of customers who are lost completely to this market, who perhaps move on to other technology or just stop using that type of product.

The simulation results are shown below for the case where all products start with 100 customers, no new customers are added, and the loss rate is zero.

Switching Based on Market Shares

There is one further variant to customer brand switching that may be useful in some circumstances (it is available by clicking here):

In this version, rather than using a switching in probability, the number of customers gained each unit time is determined by a dynamic overall market share, the one-dimensional array switching in market share.  Whereas the previous model (above) includes a mechanism whereby brand X explicitly gains customers from brand Y, this model bases switching on market share dynamics; if a customer decides to stop using brand X, they might then go searching for a replacement and would be more likely pick the market leader.  Note the switching in market share can be the same as the brand’s overall market share or, as shown, it can reflect the different market shares between new customers (market share) and switching customers (switching in market share).

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

In the session, Chris Soderquist demonstrates how you can use embeddable simulations to help raise awareness and understanding about climate change.  Using simple bathtub dynamics, the “C02 Dynamic Indicator” simulation lets anyone experiment with the current set of policy proposals for capping and reducing emissions.

The live session was recorded and is now posted on our website, along with two embedded simulations that you can run online.  All of the session materials are open source and we encourage you to link, embed and share them with others.

In the second post of this series, I showed how to selectively pull information from an array in order to route water through a watershed.  In this post, I will use the exact same technique to move customers between different product brands.

Switching Customers between Different Products

Business models often need to model gaining customers from, and losing customers to, competing products in a relatively mature market (what Kim Warren, in his excellent book Strategy Management Dynamics, calls “Type 2 Rivalry”).  These are often driven with statistical models developed through market research.  For this application, we need a matrix describing the probability of switching from product A to product B each time unit.  A sample appears in the table below.

switching probability (units: dimensionless)

To read this table, locate the product the customer is presently using in the left column (say, B).  Read across that row (the second row, in this case) until you find the product the customer is switching to (say, C).  The number in that cell (in this case, 0.05 or 5%) is the probability the customer will switch from the first product to the second (from B to C) in this time unit.  If the model is running in months, as ours is, this table indicates that 5% of customers using product B switch to product C every month.

Of course, the values in the table do not need to be constant.  Often each cell will contain a regression equation based on various product characteristics – including market share, marketing effort, product features, and product quality – that evolve over the course of the simulation.

Note the diagonal is zero.  This means customers do not switch from one product to the same product.

Note also that the sum in any row cannot exceed 1.0, which represents 100% of the customers using that product.  It is quite normal for it to be below 1.0 because we do not include people who are not switching.  Some modelers find it easier to always have each row add up to 1.0.  If you desire to do this, fill the diagonal with the difference between 1.0 and the sum of the other columns.  For example, to do this for product A, replace the top left cell with 1.0 – (0.01 + 0.03 + 0.05 + 0.001) = 0.909 [for you Beatles fans].

Exploring the Model

The model is slightly different from the watershed model in my previous post.  This model is driven by the probability of switching.  Thus, it is driven by the stock of customers for each product, while the watershed model was driven by the outflow from each reach.  This is shown below (and is available by clicking here).

The two-dimensional array switching (called a “switching matrix”) is derived by multiplying customers by the switching probability:

Customers[Products]*switching_probability[Products, Products]

This equation does not perform a traditional matrix multiplication, but multiplies each element of Customers (a one-dimensional column) by every element in the corresponding row of switching probability (two-dimensional).  [For more information about matrix arithmetic in iThink and STELLA, click here.]  As an example, given the switching probability above, if 100 customers use each of products A, B, C, D, and E, the values of switching with its row and column totals would be:

switching (units: customers)

These numbers are derived by multiplying the entire row for each product by 100, the number of customers using that product (again, this example uses 100 customers for each product).  The number of customers switching from a given brand each month is the sum of that product’s row.  The number of customers switching to a given brand each month is the sum of that product’s column.

Thus, 9.1 existing customers switch from (stop using) product A this month while 4.2 customers switch to (start using) product A this month.  Product A will therefore have a net loss of 9.1 – 4.2 = 4.9 customers for the month.

The switching matrix directly drives both the loss of customers and the gain of customers.  The former, losing_customers, sums rows:

ARRAYSUM(switching[Products, *])

while the latter, gaining_customers, sums columns:

ARRAYSUM(switching[*, Products]).

Running the Model

The simulation results are shown below for the case where all products start with 100 customers and no new customers are added.

Note the product with the largest net gain = Total In – Total Out, product D, gains the most market share over the course of the simulation, while that with the smallest net gain, product B, loses the most.

Next time, I will show a more general switching matrix formulation.

This is the second installment of a multipart series.  The first part can be found by clicking here. Part 3 is available here.

In the first post of this series, I showed how to conditionally pull information from an array.  In this post, I will extend this concept to show how to route information through an arrayed model.  This is especially useful in spatial modeling applications.

Routing Water Through a Watershed

A common ecology application is the modeling of a watershed.  Part of such a model will necessarily involve a network of stream or river segments – called reaches – which feed each other.  It is desirable to implement this in a way that makes it easy to modify the reach network.  Using an explicit stock-flow network makes this very difficult.  However, it is relatively straightforward to use arrays of stocks and flows to build an easily configurable network.

Imagine a small watershed broken down into reaches as shown below:

For our purposes, a new reach will need to be created at every junction point.  Therefore, in this example and from a topological point-of-view, it is not strictly necessary to treat reach 4 separately from reach 2 nor reach 5 separately from reach 3, but reaches 2 and 3 must be separate from reach 1.  There are, of course, other reasons to separate reach 4 from 2 and reach 5 from 3, for example, slope, channel width, length, etc.

Every reach flows into exactly one other reach at its head, but many reaches can flow into the head of the same reach.  This requires a many-to-one representation of the reach network.  This is accomplished quite easily with a routing map which, for each reach, contains the number of the reach that this reach flows into.  We also need someway to signify the outlet.  Since reach numbers start at one, we can use zero to signify the outlet.  Using these rules, the above network is completely represented in the following routing map:

The nice thing about this representation is that it fits nicely into a one-dimensional array where the array index is the reach number and the reach it flows into is the value stored in that array element.

The model itself uses one stock to represent each reach.   That stock has one inflow for water entering the reach and one outflow for water leaving the reach: (Download the zipped STELLA model here)

The Routing Magic

The outflow of one reach is directed to the next through the two-dimensional array route water.  This array is arranged with one outflow per row and one inflow per column.  Thus, the outflows are placed one per row in the proper column for the reach it is headed for.  This means there can only be one non-zero value in each row.  For example, at one point in the simulation, the array route water contains these values:

As expected, there is at most one non-zero value for each row.  There is no non-zero value in the first row because that corresponds to the outlet, so no water is routed to another reach.  The value in the second row is the outflow from reach 2, the value in the third is the outflow from reach 3, and so on.  The columns sum to the amount of water being fed to the given reach.  So reach 1 will receive 1.04 + 0.87 = 1.91, reach 2 will receive 0.33, and reach 3 will receive 0.24.  Reaches 4 and 5 receive no inflow because they are at the top of the watershed (there are no reaches feeding into them).

The equation that does this magic is:

IF (routing_map[Reaches] = ARRAYIDX(2))
THEN outflow[Reaches]
ELSE 0

This says that when the column number (ARRAYIDX(2)) equals the reach number of the downstream reach (routing_map[Reaches]), fill the cell with that reach’s outflow.  Note that we at most define one downstream reach for each reach, so this will only be non-zero in only one column.  In this way, we avoid violating conservation of water.

The inflow then needs to sum the columns to get the total flow directed into that reach:

ARRAYSUM(route_water[*,Reaches])

The outflow at the stream outlet is captured in the converter final outflow, which has the equation:

ARRAYVALUE(outflow[*], ARRAYMINIDX(routing_map[*]))

This picks up the outflow value of the reach that has the smallest entry in the routing map.  Since we defined the outlet to have the smallest value, zero, as long as we do not accidentally give two entries a zero value in the routing map, ARRAYMINIDX will find the array index (i.e., reach number) of the outlet.

With constant precipitation, the behavior of the model at the end of a storm event appears below:

Extending the Model

This particular example has two simplifications that are easy to extend.  The average time it takes water to flow through a reach (avg reach length in minutes) is constant, but can be made an array with different values for each reach.  These values can also be calculated as the simulation progresses based on existing conditions in the reach.  There is also a constant precipitation to all reaches all the time to keep water flowing into the system.  This can be replaced by an array and a more intricate model of how precipitation flows into the stream over the landscape, or by a simple time series using a graphical function.

One final note about the model:  I chose a stock with an average residence time to represent each reach.  Under some circumstances, this may generate results that are confusing.  For example, because it is an average residence time that follows an exponential distribution, some water that appears at the top of a reach will immediately exit from the bottom of the reach and other water that enters the reach will exit long after the average time has expired.  In individual circumstances, for example a single short storm event, this representation may not match the real world very well.  However, on average, over a long enough period of time, it matches the real behavior quite closely.

One alternative that enforces an exact delay for water moving across each reach is a conveyor.  However, a conveyor is discrete in its behavior, so the model will not behave as true water flow does, i.e., continuously.  Three stocks, i.e., a third-order material delay, could also be used for each reach to reduce the amount of entering water that exits immediately.  Although this does cause the model to retain more water at the start of an event, it also leads to water being lost more quickly at the end.

Next time, I will use this same structure to model customers switching between different brands.

Using arrays can be quite intimidating for most people. Many times, it is difficult to discover the correct way to formulate a problem in terms of arrays, especially when trying to do so in terms of single equations that can be applied to all elements of the array.

Consider the case where you might wish to count the number of occurrences of a value in an array. This can arise in many applications that need to track attributes, but is prevalent in spatially-explicit business applications. In such an application, you may associate a product code with a location and then want a count of the products of a given type. The following examples demonstrate a common way to extract conditional information from an array.

Finding the Number of Stations with a Given Status

Imagine you have a two-dimensional grid of fire stations in a city, called Stations, that stores one of four statuses:

0: no station in this sector
1: ready
2: away on a call
3: refitting

You consider the number of fire stations ready at any given moment to be an important metric. To calculate this, connect Stations to another two-dimensional array of the same size called ready stations. This will have a one in an array element if the station for that quadrant is ready and a zero otherwise. Its equation is:

IF Stations[Y, X] = 1 THEN 1 ELSE 0 { station ready? }

Note this equations uses dimension names (i.e., Stations[Y, X]) rather than element names (e.g., Stations[1, 2]). This allows you to create just one equation for the entire array (with “Apply to All” turned on), rather than a separate equation for each individual array element (with “Apply to All” turned off). When “Apply to All” is turned on, the equation for each element of the array is automatically generated by substituting that element’s dimensions for the dimension name in the given equation. All of the examples in this post use dimension names.

The total number of ready stations is now just the sum of all of the elements in the array ready stations. This is easily calculated by connecting ready stations to a scalar converter named total ready stations that has the equation

ARRAYSUM(ready_stations[*, *])

The model is shown below and can be downloaded by clicking here.

This general method can be used anytime we need to count the number of elements in an array based on some condition. First, create an array that has its elements set to one if the array-based condition is met and zero otherwise (IF condition THEN 1 ELSE 0). Then create a converter to sum the elements of this new array (using ARRAYSUM). Remember to turn “Apply to All” on and use dimension names in the condition rather than element names.

Finding the Number of Stations with Each Status

We can extend this method to map data from an array indexed by one attribute (or dimension), say location, to another array indexed by a different attribute, say status. Imagine now that instead of just counting how many stations are ready, you want to know how many stations have each status.

One way to do this is to extend the structure above to separately count the number of stations in each status, as shown below.

In this case, away stations and refitting stations have the following formulas, respectively:

IF Stations[Y, X] = 2 THEN 1 ELSE 0 { station away? }
IF Stations[Y, X] = 3 THEN 1 ELSE 0 { station refitting? }

Walking Through a Specific Example

Let’s look at what happens in this model if array of fire stations is only one-dimensional, Stations[Location], and each location is named. This eliminates the “no station in this sector” status. Consider a five-location fire department with the following statuses in each location, held in the array Stations:

The structure above would generate the following values in the three separate arrays ready stations, away stations, and refitting stations. Their respective sums, in the three corresponding “total status stations” converters, are shown at the bottom.

While this method works, it does not have a very flexible structure. It requires structural changes whenever a new status is added (or removed) and also needs to be customized for every individual application. There are also applications (for example, tracking products) where it is necessary to move between arrays that are indexed by location, as Stations is, and arrays that are indexed by the data within the location, in this case by status.

A More General Solution

If we define another array dimension called Status, with element names Ready, Away, and Refitting that correspond to the numeric statuses 1, 2, and 3, we can represent this same solution in one two-dimensional array station map[Location, Status]:

The array station map allows us to map stations directly to their statuses. Each row of the map represents one station. By definition, there can only be one “1″ in each row because a station can only have one status. The sum of each row is the number of statuses that are assigned to that station, which must always be one. Each column of the map represents one status. Multiple “1″s can appear in the columns because many stations can have the same status. The sum of each column is the number of stations assigned to that status, the value we are trying to determine.

Note the similarity between this two-dimensional array and the preceding table that shows the contents of the individual one-dimensional arrays ready stations, away stations, and refitting stations. Each column of station map corresponds to one of these one-dimensional arrays.

The model for this more general solution is shown below and can be downloaded by clicking here.

The equation for station map[Location, Status], which will map Stations as illustrated above, is

IF Stations[Location] = ARRAYIDX(2) THEN 1 ELSE 0

Note the use of ARRAYIDX() to parameterize the status we are interested in for this array element. The parameter to ARRAYIDX, 2, gives the index of the second dimension for this array element, i.e., the index that corresponds to Status. Since Apply to All is on, this will become the desired status for each array element.

The number of stations with each status is calculated in the one-dimensional converter station totals[Status] using the equation

ARRAYSUM(station_map[*, Status])

This equation sums the columns of station map. For the above example, station totals will be filled with the values 3 (Ready), 1 (Away), and 1 (Refitting). Note we have used an intermediate two-dimensional array, station map, to properly map our data from a one-dimensional array indexed by the attribute Location to another one-dimensional array indexed by the attribute Status.

Next time, I will expand this concept to create general routing structures.

• Navigation buttons are styled to look like hyperlinks instead of push buttons.
• Push buttons are used to perform actions, such as running the simulation.
• Tabs along the top of each page let the viewer know which section they are in, as well as providing a way to quickly jump to any section.

An example of the interface template applied to a model published online

You can download the template model and use it as a starting point for creating interfaces to models you want to publish online.

Download the template for STELLA

Download the template for iThink

How to use the template

The template model consists of four interface pages: Home, Background, Explore Model and Simulate.  Navigation buttons to and from each page are already setup within the interface.  You can customize each section by replacing the placeholder text.  The pages are a starting point that establishes a website-like look and feel.

The easiest way to apply the template model’s interface to your own model is to simply copy and paste the pages.

1. Open your own model and switch to the Interface layer.  Choose ‘Interface Prefs’ from the Interface menu and make sure the ‘Show Pages’ option is selected.  Set Page Size to 800×600.  This will allow you to see the page boundaries as you apply the template.
2. Open the template model and from the Interface layer choose the Edit menu –> Select All, then Edit menu –> Copy to copy all of the pages to the clipboard.
3. Go back to the interface layer of your model and choose Edit menu –> Paste.

When you have pasted the pages into your model, they will probably not be lined up correctly within the page boundaries.  The dashed lines show the boundaries of each page.  Using these dashed lines as a guide, use the arrow keys on your keyboard to make the pasted pages line up properly.

Now you can customize the title and text in each page.  Add a picture to the home page and build a control panel.  If you want, you can use the hyperlink button on the ‘Background’ page to link to online information, (PDFs, other models, Wikipedia, etc.).  Consider building a storytelling sequence on the ‘Explore Model’ page to unfurl the model a chunk at a time.  These are just a few customizations you can make.  Take a look at some of our sample models published online with isee NetSim for ideas.

If you publish some models online that use the template, we’d love to see them!  Leave a link in the comments so we can check them out!

Whether you think of the program as a stimulus for the economy or the environment, it cannot continue to operate without funds.  Congress had to quickly approve another $2 billion for the program when the initial $1 billion appropriation was exhausted in less than a week.  The rebate fund lacks an inflow of cash to offset the outflow of payments.

I wrote about a STELLA model of a self-financing program to promote cleaner vehicle sales in May.  The “Feebates” model uses fees charged on sales of vehicles that have high pollution rates to finance rebates for cleaner vehicles.  I asked the Feebates model author, Andrew Ford what he thought about “Cash for Clunkers”.  He sent me an article he co-authored with Todd BenDor titled “Simulating a combination of feebates and scrappage incentives to reduce automobile emissions”.  Basically, the article describes a feebates model with a “cash for clunkers” type of rebate added to promote scrapping older vehicles.  The article was published in the journal Energy, 31(2006), 1197-1214.

Mr. BenDor was kind enough to send me the STELLA model so we could publish it to the web with isee Netsim.  We thought this would be a great way to start some discussion around the Clunkers program, public policy and using simulations to facilitate decision making.

Click here to run the simulation online.

• Can you design a policy that stimulates the economy and reduces emissions?
• This simulation uses fees to fund the program instead of deficit spending.  What other sources of revenue could we consider?
• Many states levy a “sin tax” on tobacco and alcohol.  How would the public react to fees on high-pollution vehicles?

Post your answers in the comment section below.

Image credit: http://www.flickr.com/photos/threadedthoughts/3808715790/

Marcus Robinson & James Ranney accept 2009 Barry Richmond Scholarship Award

The Barry Richmond Scholarship Award was established in 2007 by isee systems to honor and continue the legacy of our founder, Barry Richmond.  Barry was devoted to helping others become better “Systems Citizens”.  It was his mission to make Systems Thinking and System Dynamics accessible to people of all ages and in all fields.

The award is presented annually at the International System Dynamics Conference to an individual whose work demonstrates a desire to expand the field or to apply it to current social issues.

Barry had a particular passion for expanding the field into K-12 education.  He worked closely with teachers and young people at all levels to incorporate Systems Thinking tools and methodologies into the classroom.  With this in mind, it was a great pleasure to present this year’s award to Marcus Robinson and James Ranney, two students from Wilson High School in Portland, Oregon.

Marcus Robinson presents at SD Conference

The scholarship committee was very impressed by the quality of the models and technical papers that these students submitted.   Both high school seniors were enrolled in Diana Fisher’s modeling class this past year.  For their final project assignment, they were asked to explore a topic of their own choosing.

Inspired by a video of the PBS NOVA program “Dimming the Sun”, Marcus chose to create a model that would help to answer the question “How Much do Carbon Emissions Need to be Reduced in Order to Stop Global Warming?” To download Marcus’ paper, click here.

After a run-in with a large frog at his summer camp, James was interested in answering “What is the Impact of Introducing an Invasive Species into an Ecosystem?” – in this case, the bullfrog into the Pacific Northwest. To download James’ paper, click here.

James Ranney presents during K-12 poster session

Both students presented their work at a parallel session during the conference and again at a special K-12 poster session.  What a delight it was to see how empowering modeling real-world issues is for young people and the enthusiasm they share for their work.

Thank you to the society for making a special effort to serve the K-12 community and provide opportunities for partnership.  Congratulations Marcus and James, Barry would be very proud!

Applications for the 2010 Barry Richmond Scholarship Award will be available on the isee systems and System Dynamics Society web sites. Check those sites for more information.

We hope to see you at as many of the following isee sponsored events as your schedule allows:

Monday, July 27th

Tuesday, July 28th

Wednesday, July 29th

Thursday, July 30th

See you in Albuquerque!

What I’ve come to appreciate since then is how the characteristic themes of the Systems Archetypes transform across all sorts of different fields and situations — even our personal lives!

Take, for example, the basic story line of the “Success to the Successful” archetype:

When given the choice, we invest our resources where we expect them to deliver the best results.  By giving more resources to one option over another, we create a self-fulfilling prophesy whereby the favored option is perpetually more successful.

The story of the Success to the Successful archetype applies to all sorts of  situations leading to well-known patterns of behavior:

• rich getting richer and poor getting poorer
• winners being rewarded with the means to win again
• popular brands attracting more popularity
• one species driving another to extinction
• top salespeople being paid higher commissions

Exploring the Model Structure

We recently published a model of the Success to the Successful archetype to the web using isee NetSim. Exploring the model is a great way to understand the underlying structure of the system and think about ways to avoid the problems it creates.

You’ll also get an appreciation of how the decision policy for allocating resources can determine success rather than competence.

Running the Simulation

After you’ve explored the model, try running a simulation.  The base case scenario assumes no one has an advantage over the other.  As you can imagine, everyone is equally successful and it’s a win-win situation.  Try turning on the “Advantage A Switch”  to see how even a modest advantage for A can snowball into a disadvantage for B.  It’s surprising how quickly the gap can widen.

Using Modules to Create Models

In STELLA and iThink version 9.1, we added the ability to build models by linking together modules.  The Success to the Successful model is an example of how you can use modules to create a higher level map of your model.  This map can easily be presented as a causal loop diagram.

The beauty of modules is they simplify the process of transitioning from a CLD to a model that actually simulates.  If you’ve ever tried to convert a causal loop diagram into a stock and flow model, you can appreciate what I’m talking about!  By architecting your model into modules, you’ve got a built-in mechanism for developing your model in manageable chunks and communicating the high level causal relationships.

• Addition and Subtraction
• Transposition
• Extracting the Diagonal of a Square Matrix
• Multiplying and Dividing by a Scalar
• Multiplying a Square Matrix by a Vector
• Multiplying a Non-Square Matrix by a Vector
• Multiplying Two Matrices

In the examples provided, I use five different arrays of varying dimensions. They are all defined using the subscript name Dim3, which is size 3.

a: 3×3
b: 3×3
c: 3×3
d: 3×1 [one-dimensional]
e: 3×1 [one-dimensional]

I also use two arrays based on Dim3, as well as the subscript name Dim2, which is size 2.

s: 2×3
t: 3×2

Addition and Subtraction

The easiest functions to perform are addition and subtraction. When two arrays have the same dimensions, they can be added and subtracted in the usual way, placing the result in a third array of the same dimensions.

For example, to place the difference of b and c into array a, set the equation of a to

b[Dim3, Dim3] – c[Dim3, Dim3]

Transposition

Transposition of a 2-dimensional array is also relatively easy, especially when the dimensions differ. Suppose we wish to transpose s and place the result in t. Simply set the equation for t to

s[Dim2, Dim3]

This will properly fill t, which has dimensions Dim3xDim2, with:

s₁₁   s₂₁
s₁₂   s₂₂
s₁₃   s₂₃

However, if the two dimensions of the array have the same name, the ARRAYVALUE built-in must be used. Suppose we wish to transpose b, placing its transposition in a. Set the equation of a to

ARRAYVALUE(b[*, *], ARRAYIDX(2), ARRAYIDX())

Note the index of a’s second dimension is used as the index for b’s first dimension and the index of a’s first dimension is used as the index for b’s second dimension.

Extracting the Diagonal of a Square Matrix

Extracting the diagonal from a 2-dimensional array is also very easy. To place the diagonal of a into d, set the equation of d to

a[Dim3, Dim3]

Multiplying and Dividing by a Scalar

Multiplying or dividing by a scalar works as expected. The contents of b can be multiplied by a scalar x, for example, by setting the equation of a to

x*b[Dim3, Dim3]

Multiplying a Square Matrix by a Vector

Matrix multiplication, as defined in linear algebra, is not directly supported at present. It is possible to multiply two arrays, but the results will not match matrix multiplication. When you multiply a 1-dimensional array by a 2-dimensional array, the size of the 1-dimensional array must match either the number of rows or the number of columns in the 2-dimensional array. In both cases, iThink and STELLA treats the 1-dimensional array as either a column or row of scalars that are multiplied across each corresponding row or column, respectively, of the 2-dimensional array.

Given the following equation in a,

d[Dim3]*b[Dim3, Dim3]

d will be treated as a column of scalars (when the 2-dimensional array is square with the same dimension name in both directions, the 1-dimensional array is always treated as a column vector). Each scalar in d will be multiplied across each corresponding row of b. Note this equation does not represent a valid matrix multiplication, and also that switching the order to b[Dim3, Dim3]*d[Dim3], while now a valid matrix multiplication, will not change the result. I recommend leaving the order as above to emphasize that it is not a normal matrix multiplication.

More precisely, the array a will be populated with these values:

d₁*b₁₁   d₁*b₁₂   d₁*b₁₃
d₂*b₂₁   d₂*b₂₂   d₂*b₂₃
d₃*b₃₁   d₃*b₃₂   d₃*b₃₃

Note that these are not the correct terms for matrix multiplication. If we multiply b[Dim3, Dim3]*d[Dim3], which is the only correct way to express it, the linear algebra result is a one-dimensional vector with these values:

d₁*b₁₁ + d₂*b₁₂ + d₃*b₁₃
d₁*b₂₁ + d₂*b₂₂ + d₃*b₂₃
d₁*b₃₁ + d₂*b₃₂ + d₃*b₃₃

If we desire to generate these terms, we must first transpose b by placing the transposition equation given above in c:

ARRAYVALUE(b[*, *], ARRAYIDX(2), ARRAYIDX())

This puts the following values in c, i.e., b transposed:

b₁₁   b₂₁   b₃₁
b₁₂   b₂₂   b₃₂
b₁₃   b₂₃   b₃₃

The equation for a can now be set to

d[Dim3]*c[Dim3, Dim3]

This equation will populate a with these values:

d₁*b₁₁   d₁*b₂₁   d₁*b₃₁
d₂*b₁₂   d₂*b₂₂   d₂*b₃₂
d₃*b₁₃   d₃*b₂₃    d₃*b₃₃

Note that this gives the proper terms for the matrix multiplication, arranged vertically. The final answer is found by summing each column, using the following equation in vector e:

ARRAYSUM(a[*, Dim3])

Aside: The terms we get from the original equation for a, d[Dim3]*b[Dim3, Dim3], are the correct terms for d^Tb, i.e., d-transposed times b, again arranged vertically. If this was the desired result, we do not need to first transpose b and can just sum the columns into vector e as shown above.

Multiplying a Non-Square Matrix by a Vector

If the matrices do not use the same dimension name for both dimensions (square or not), you have much more control over the resultant terms and do not have to resort to transposition at all. The following discussion uses the following additional arrays:

f: 2×3
g: 2×1 [one-dimensional]
s: 2×3t: 3×2

If we multiply f times d, we will get a different result based on whether the resultant array is 2×3 or 3×2. Placing the equation d[Dim3]*f [Dim2, Dim3] into array s yields exactly the terms needed for a matrix multiplication, arranged horizontally:

d₁*f₁₁   d₂*f₁₂   d₃*f₁₃
d₁*f₂₁   d₂*f₂₂   d₃*f₂₃

Setting the equation of g to

ARRAYSUM(s[Dim2, *])

will properly sum the rows into g, giving the result of the matrix multiplication f *d.

If, instead, we place the same equation d[Dim3]*f [Dim2, Dim3] into array t, we get these terms instead:

d₁*f₁₁   d₁*f₂₁
d₂*f₁₂   d₂*f₂₂
d₃*f₁₃   d₃*f₂₃

The correct terms are still available, only now they are arranged vertically and the columns must be summed to arrive at the proper result:

ARRAYSUM(t[*, Dim2])

Multiplying by g instead, i.e., placing the equation g[Dim2]*f [Dim2, Dim3] into s, yields these terms:

g₁*f₁₁   g₁*f₁₂   g₁*f₁₃
g₂*f₂₁   g₂*f₂₂   g₂*f₂₃

These are the terms for g^Tf , arranged vertically. Placing the equation into t gives the same terms horizontally. To calculate f ^Tg, it is necessary to first transpose f by setting t to the equation

f [Dim2, Dim3]

When this is multiplied by g, the correct terms for f ^Tg will be generated (horizontally in a 3×2 array and vertically in a 2×3 array).

Multiplying Two Matrices

Unfortunately, there is no way to easily generate the terms for multiplying two matrices. The equations must be individually tailored for each case. As an example, let’s look at multiplying t*s and storing the result in a. Array t is 3×2 and array s is 2×3, so the result will fit in a (3×3). The correct answer is:

t₁₁*s₁₁ + t₁₂*s₂₁    t₁₁*s₁₂ + t₁₂*s₂₂    t₁₁*s₁₃ + t₁₂*s₂₃
t₂₁*s₁₁ + t₂₂*s₂₁    t₂₁*s₁₂ + t₂₂*s₂₂    t₂₁*s₁₃ + t₂₂*s₂₃
t₃₁*s₁₁ + t₃₂*s₂₁    t₃₁*s₁₂ + t₃₂*s₂₂    t₃₁*s₁₃ + t₃₂*s₂₃

Examining these values for a pattern that could be used in an apply-to-all equation, we can see the proper equation for array a is

t[Dim3, 1]*s[1, Dim3] + t[Dim3, 2]*s[2, Dim3]

This pattern, incidentally, will hold in general for any two matrices. However, this equation does not work as an apply-to-all equation because the same name is used for both the first and second dimensions so it is not possible to tell which of the two varying Dim3 values to use. In this case, iThink and STELLA currently always substitutes the value of the first dimension, which will not give the correct answer. This will always be a problem when the resulting array is square. One workaround is to define an additional dimension name that is the same size as Dim3. This works, but may lead to other array incompatibilities. The other solution is similar to how we transposed a square matrix – explicitly use ARRAYVALUE for any references to the second dimension:

t[Dim3, 1]*ARRAYVALUE(s[*, *], 1, ARRAYIDX(2))
+ t[Dim3, 2]*ARRAYVALUE(s[*, *], 2, ARRAYIDX(2))

This is certainly not as clear, but it works and is only necessary when the result is square. One last example underscores this point and also shows how to expand the general pattern to a larger array. For this example, one more dimension name, Dim4, of size 4 and two more arrays are needed:

u: 3×4
v: 2×4

To calculate the result of s*u in array v, where s is 2×3, u is 3×4, and v is (correctly) 2×4, set the equation for v to:

s[Dim2, 1]*u[1, Dim4] + s[Dim2, 2]*u[2, Dim4] + s[Dim2, 3]*u[3, Dim4]

Note this same method can be used for multiplying a 2-dimensional array by a 1-dimensional array instead of the method presented above. It avoids both the intermediate transposition and product arrays, but results in more complex expressions that must be tailored to each combination of array sizes.