Using system dynamics for long term bottom-up electric load modeling in
rural electrification
Elias K. Hartvigsson Jimmy Ehnberg
Electric Power Engineering, Electric Power Engineering,
Elias.hartvigsson@chalmers.se
Erik O. Ahlgren Sverker Molander
Energy Technology Environmental Systems Analysis
Chalmers University of Technology
Abstract: One billion people lack access to electricity, most of whom live in rural areas in
developing countries. A solution to improve their situation is small independent grids, so
called minigrids. One of the major challenges for minigrids is to become economically
reliable. As electricity usage is a major source of income for a utility, it is important to
consider how its fluctuations impacts a utility. This work presents an integration of a
previously developed system dynamics model with a comparably detailed bottom-up load
model developed in MATLAB. The simulations show that while using a more detailed load
model results in an increase in generation capacity expansion frequency and that the
investments are made in smaller sizes. Due to the different approach to integrating
electricity usage growth, the bottom-up load model shows a faster increase in electricity
usage than the system dynamics load model. With a modeled difference on net income,
power utilization rate and electricity usage the results indicate that it is important to
consider improved load model detail when modeling income and expenditures of an electric
utility.
Keywords: bottom-up load modeling, system dynamics, mingrids, rural electrification
Introduction
As of today one billion people lack access to electricity around the world. Roughly half of
these people live in sub-Saharan Africa, and a large majority of them live in rural inaccessible
areas (IEA, 2013). Improving electricity access for these people is considered an important
step towards combating poverty and improve their social and economical conditions. As a
service, access to electricity is considered to be a necessary but not sufficient condition for
economic and social development (Barnes, 2007; Bhattacharyya, 2013; Goldemberg,
Johansson, Reddy, & Williams, 1985). In many cases, using traditional methods such as grid-
extension is not enough and if rural communities are to gain the full benefits of electricity
access within the foreseeable future, reliable and well managed off-grid solutions are
needed. (Ahlborg & Hammar, 2014; Diaz, Arias, Pefia, & Sandoval, 2010; IEA, 2013;
Tenenbaum, Greacen, Siyambalapitya, & Knuckles, 2014; Urpelainen, 2014)
The size of off-grid solutions varies from tens of W for small Solar Home Systems (SHS) to
large minigrids with generation capacity of several MW. The low costs of SHS makes them
affordable for single households but their capacity is limited to supply a few low consuming
appliances such as: lights and radios. Thereby limiting the benefits that can be obtained from
SHS (Azimoh, Klintenberg, Wallin, Karlsson, & Mbohwa, 2016). Minigrids can on the other
hand supply enough power and energy to be used for productive activities such as milling,
workshops, irrigation pumps, welding and shops/bars. But needs a managing organization
and the funneling of resources. Minigrids have been used in rural electrification with various
levels of success, with one of the major challenges being their poor economic performance,
leading to an inability to reach cost-recovery (Barnes & Foley, 2004; Kirubi, Jacobson,
Kammen, & Mills, 2009; Levin & Thomas, 2014; Schnitzer et al., 2014). Their ability to reach
cost-recovery has been related to the systems utilization factor (Kirubi et al., 2009; Sarangi
et al., 2014). A systems utilization factor is the fraction of actual to maximum electricity
generation. In order to improve the utilization factor, generation should be matched with
current demand and be appropriately adapted to future demands.
Electric load modeling is a common tool in order to analyze current and future electricity
demand. It can be divided into two types: top-down and bottom-up modeling. Top-down
load modeling is based on large aggregated data sets and use mathematical methods to
analyze and predict demand changes for large groups. Bottom-up load modeling, on the
other hand, is based on modeling of the single loads, such as lights and TVs, which are then
aggregated to obtain a total user and system load.
One advantage of top-down methods is that they require less consumer data than bottom-
up methods. Bucher and Andersson (2012) developed a top-down approach to load profile
generation for household load profiles. By using a top-down approach their method could be
applied in cases where little or no consumer information is available. However, in general,
top-down approaches used for residential loads are less frequent than bottom up methods.
Alternatively, they are more often used to aggregated load analyses since their accuracy
decrease further down in the system, partly due to lack of integration with human behavior
and specific appliance usage.
Bottom-up methods have the benefit that they either directly or indirectly incorporate
human behavior (Stokes, 2005). This can be done by defining activities, and then generating
load profiles from these activities (Widén et al., 2009). In the case of Widén et al, they
addressed e.g. vacuum cleaning but in developing countries it could include activities such as
milling, welding or studying. The challenge with bottom-up load modeling methods is that
they are relying on large amounts of information about each specific load; data which is
often protected by integrity issues or by privacy laws.
Regardless of the choice of method, if the model simulates current electricity consumption
they do not likely need to take into account changing socio-economic conditions like
electricity price, access to new appliances, income or population. However, if the purpose is
to model future consumption with a sufficiently long time scale, local conditions change
which needs to be incorporated into the modeling process. Recently, artificial neural
networks (Carpinteiro, Leme, de Souza, Pinheiro, & Moreira, 2007), particle-swarm-
optimization (Niu, Li, Li, & Liu, 2009) and other artificial intelligence methods have been
employed to predict the long term impact of changing conditions on load changes in the
western world.
In rural communities in developing countries, where electricity has recently been
introduced, it is likely to believe that the time scales for local conditions to change are
shorter. Introduction of new electrical appliances can have a substantial impact on electricity
provision and livelihoods. E.g. low energy lightbulbs that can reduce electricity consumption
and introduction of fridges and freezers, which both have a large impact on electricity
consumption patterns and on life style (Madubansi & Shackleton, 2006). Therefore, it is
important to consider changing local conditions in electricity load modeling in rural
electrification.
System dynamics has been explained as a tool to map system structures to system behavior
(Davidsen, 1992). As such it has been used to describe problems related to aggregated socio-
economic indicators, such as income, energy usage and environmental impact (Hjorth &
Bagheri, 2006; Marzouk & Azab, 2014; Qudrat-Ullah & Davidsen, 2001). In the case of energy
modeling in general and electricity modelling in particular, there have been several
applications of system dynamics, see Teufel et al. (2013) for a review. However, modeling
the elements of aggregated variables falls outside the scope of system dynamics. In the case
of electricity usage, system dynamics can be an efficient tool to analyze growth in
aggregated electricity consumption but is not suited to model daily load profiles.
A system dynamics model of a utility operating a rural minigrid in a developed country
(Tanzania) was presented in (Hartvigsson, 2016). The model focus on the income and
expenses of a utility, which is directly linked to electricity usage and generation capacity.
However, the model assumes a fixed, linear relationship between electricity consumption
and peak load. In reality, the fast daily variations of different user groups are much more
complex. Inclusion of these fast variations is likely to result in higher peak loads and,
thereby, a higher peak power demand for a given electricity usage. To the knowledge of the
authors, no link between a system dynamics electric utility load model and a load model has
previously been presented in the literature. Despite the self-evident difference between the
two load modeling approaches, the impact and its extent is not known. Therefore, the
purpose of this paper is to investigate to which extent improved load modeling impacts a
system dynamics utility model. Specifically, does improved load modeling detail impacts a
modeled utility’s economic performance? Economic performance is defined as the utility’s
weekly net income during the simulated time period, e.g. differences between income and
expenses.
The paper is divided into five sections. First the method of integrating the load model with
the system dynamics model is presented. Second, the load model is explained in detail
together with input data and assumptions. Based on the two model approaches, a number
of runs are shown both with and without an external load model in the Result section. The
results and their possible impacts are discussed before conclusions and future work are
presented.
Method
This work is based on the system dynamics model presented in (Hartvigsson, 2016). The
model employs an aggregated and simplified approach to load modeling (from here on
referred to as the Vensim load model), where peak demand is scaled from electricity usage
electricity usage
da
according to peak demand = -, where d is the number of hours of electricity
consumption per week. In order to investigate the impact of increased load modeling detail
has on the modelled utility’s economy, a stand-alone load model was developed in MATLAB
(from here on referred to as the MATLAB load model). The results on the modelled utility’s
economy and overall electricity usage are then compared between the simplified system
dynamics load model and the MATLAB stand-alone load model.
The system dynamics model was developed in Vensim, making it possible to utilize the
VenDLL library for handling external calls. The VenDLL library allows external calls to change
selected variables during a Vensim gaming simulation. This allows VenDLL to retrieve data
from the Vensim simulation, send the data to the load model in MATLAB, run the load model
and retrieve its output and then write the changes to Vensim. The communication between
Vensim and MATLAB is conducted at each time step, allowing for real time integration of the
MATLAB load model into the system dynamics model. Figure 1 shows a conceptual diagram
of the modeling process.
- Power Availability
- Electricity Usage
Load Profile
MATLAB
function [y] = scaleHHTV(x,D)
k=0.002;
- Income
- Number of users
Figure 1 A conceptual diagram of the integration of Vensim with the load model in MATLAB.
The system dynamics model consists of five sectors: user diffusion, utility economics, local
market, population and electricity generation. By including the different sectors, the model
can give a simplified representation of local socio-economic conditions that directly and
indirectly affect electricity usage. It is assumed that the main driver for increase in electricity
consumption is income. Income is sent from the system dynamics model in Vensim to the
MATLAB load model, which use it to generate individual user load profiles (see the below
load model description). The user’s load profiles are then scaled with the number of users to
obtain the total system load profile.
Both load modeling approaches assumes that if there is an initial generation capacity
installed and that the this can continuously produce at its rated power. Depending on the
energy source this might not always hold true. Specifically, renewable energy sources are
often characterized by having a relative large intermittency. However, modeling individual
energy sources characteristics is outside the scope of this work.
Load Model
The implemented load model is a linear bottom-up appliance diffusion model using income
and power utilization rate as variables. As a bottom-up load model it models each appliance
separately and then aggregates them to generate the system’s total load profile. In order to
determine the load growth for each appliance, they are assigned two variables: a diffusion
rate and a saturation level. The diffusion rate determines how sensitive appliance load
growth is to variations in income and the saturation limit marks the maximum occurrence of
each appliance. Equation 1 describe the load profile for appliance i.
I, =1+8+StdLP;-N(u,9) (1)
Where | is average income, 6 is diffusion rate, StdLP; is the standard load profile for
appliance / and N(u, ~) is a normal distribution with mean py and standard deviation g. By
selecting mean and standard deviation to be 1 and 0.1, the model allows for some
uncertainty of each appliance power output at a specific time interval. This reflects part of
the uncertainty in when a specific appliance is turned on. In order to obtain the systems
total load, the appliance load profiles are aggregated first per user and then for all users.
Equation 2 then describes a generic expression for the total load profile.
Wh DEa Up T; * f Power_Availbility) (2)
Where m is the number of user groups, n the number of appliances in each user group and
U, the number of users in user group k. Due to differences in consumption patterns,
sensitivity to power outages and power demand, users where separated into two groups:
households and Small and Medium Enterprises (SMEs). Users are also assumed to respond
negatively to disruptions in the electricity provision, which is expressed by the function f
with power utilization rate as only variable. Power utilization rate is defined as the ratio
between peak power and maximum generation capacity. If this ratio becomes larger than
one, implying that temporary blackouts occur, users are assumed to respond by reducing
their electricity consumption. The function f is expressed as a cotangents function, assuming
a slow initial reaction that increase as power utilization rate deteriorates.
Load model data and assumptions
The standard load profile for each appliance was identified via a case study in a minigrid in
Tanzania (Hartvigsson, Ehnberg, Ahigren, & Molander, 2015) and rely on both interviews and
measurements. The interviews identified four appliances in households: lights, stereo, TVs
and DVDs. Some other appliances such as computers and iron where also identified but
where concluded to have no impact on average daily load profiles. Figure 2 shows the
standard load profiles for each of the appliances using a 30 min time resolution.
Lights
150 30 Stereo
= =
£ =
= 220
2 ra
g g
3 B10
3 3
i 2
we Ww
0
12 18 0 12 18 24
Hour of the day Hour of the day
Tv DvD
6 4
= =
= =
3
Ba =
g Be
2 3
5 51
2 2
ai a
0 0
0 6 12 18 24 0 18 24
Hour of the day
12
Hour of the day
Figure 2 showing standard load profiles for each of the household appliances.
For SMEs a wider range of appliances where identified. Due to their similarity in terms of
average electricity consumption some of them where grouped together. This includes the
groups “Machines” and “Others” seen below in Figure 3. Machines is a combined load group
of high power electric equipment such as electric machines, power tools and electric welding
equipment. The group “Others” mainly include low load appliances (such as hair trimmers)
or appliances that are very rare (computers and printers). In Figure 3 shows the standard
SME load profiles for each appliance.
150 Lights Stereo
25
z 20
100 =
= =45
8 8
5 50 30
2 25
Ww w
0 0
0 12 18 0 12 18
Hour of the day Hour of the day
6 Machines 6 Other
= =
£300 15
= =
° 2
8 200 B10
> >
3 3
$ 100 $5
2 2
a a
0 0
0 24 0 24
12 18 12 18
Hour of the day Hour of the day
Figure 3 shows standard load profile for each of the SMEs appliances.
As shown in Equation 1, the load profile expression is a linear function of average income.
Assuming that no income results in no appliance ownership one data point is required for
determine diffusion rate. The diffusion rate for each appliance could therefore be
determined using data on income of 29 households and 19 SMEs from a minigrid in rural
Tanzania.
The saturation level, which is defined as the maximum occurrence of a specific appliance,
e.g. a saturation limit of 14 implies 14 is the maximum amount of appliances per user
regardless of income. In some cases data on saturation levels where found. In cases where
data was lacking, values were estimated based on observations on appliance occurence and
income levels. Diffusion rates and saturation levels for all appliances can be seen in table 1.
Table 1 Showing diffusion rates and ion limits for for each applian
Diffusion rate Saturation Limit
Appliances [appliances/USD] [# appliances]
Lights (household) 0.39 14*
Stereo (households) 0.02 2
TV (households) 0.004 2
DVD (households) 0.03 1
Lights (SME) 0.03 4
Stereo (SME) 0.013 1
Electric Machines (SME) 0.002 0.5
Other (SME) 0.008 1
" Data taken on saturation limits for lights in Argentina and Sri Lanka (McNeil, 2008).
? Estimated using data from (McNeil, 2008).
Results
The results of running the Vensim and MATLAB load models are presented in Figure 4-7 in
the following order. Figure 4 shows the MATLAB and Vensim generated load profiles at
initial and final stages. Figure 5 presents the utility’s balance sheet and Figure 6 the power
utilization rate. Finally, total electricity usage is shown in Figure 7.
Figure 4 shows initial and final load profiles. The dark and light orange load profiles are
generated using the MATLAB load model as described in this paper. The dark and light blue
lines are generated from Vensim’s peak load model, as described in the Method section.
Since the load implementation in Vensim is a constant scale between electricity usage and
peak load, the Vensim generated peak loads are constant. The solid dark blue line represents
the initial Vensim load and the dotted light blue line represents the final Vensim load. The
orange load profile was generated by the MATLAB load model. The dark orange graph is the
initial load profile and the light orange graph is the final load profile. As seen there is a
distinct growth in peak power and electricity usage for the two models, but with a
considerable higher peak power shown by the MATLAB load model.
450 1
inal load profile - Matlab load model
itial load profile - Matlab load model
400} | peak load Vensim load model 4
initial peak load Vensim load model
234 5 6 8 9 10 11 12 13 14 15
Hour of the day
Figure 4 The graph shows initial and final load profiles using the external load model developed in MATLAB and the internal
Vensim load model.
Figure 5 shows the utility’s net income plotted as a function of time. The blue graph, which
represents the Vensim load model shows fewer but larger dips. The dips are explained by
sudden large expenses, and is in this case related to the acquisition of new generation
capacity. As investments in new capacity depend on power utilization rate, the expenses
follow power utilization rate with a delay. As seen in the case with the MATLAB load model,
there are more but smaller dips, indicating that expansion is taken place more often but and
in smaller sizes.
450 T T T T T T T T
fensim load model
400} jatlab load model 4
350}
8
8
Net income [k$]
8 8
6 6
a
3
100/
oO 100 200 300 40 600 700 800 900 1000
0 500
Time [Week]
Figure 5 The graph shows the utility's financial balance using the external load model developed in MATLAB (blue graph)
and the internal Vensim load model (orange).
Figure 6 shows the power utilization rate plotted over time for the Vensim load model (blue
graph) and the MATLAB load model (orange graph). The figure clearly shows a higher overall
power utilization rate when using the MATLAB load model. The MATLAB load model is also
seen to be showing less fluctuations than the Vensim load model.
1.5
fensim loa el
Matlab load model
°
Power utilization rate [p.u.]
oO 100 200 300 400 500 600 700 800 900 1000
Time [Week]
Figure 6 The graph shows power utilization rate, defined as fraction between peak power demand and installed generation
capacity, for the external load model developed in MATLAB (orange graph) and the internal Vensim load model (blue graph).
Figure 7 shows total electricity usage for the two models. The blue graph shows total
electricity usage using the Vensim load model and the orange graph using the MATLAB load
model. The total electricity usage is very similar initially and it is not until after about 300
weeks that they start to diverge with a final difference of about 25%. Even though number of
users are not presented in detail in the paper, the number of users are the same for both
cases, indicating that it is electricity usage per person that grows faster when using the
MATLAB load model.
35 r T T r " u i
Vensim load model
Matlab load model
30+
3
25
zs
=
=
=20
°
&
8
215
2
S
B10
2
Ww
Se 4
0 1 1
oO 100 200 300 400 600 700 800 900 1000
500
Time [Week]
Figure 7 The graph shows total electricity usage in the minigrid for the external load model developed in MATLAB (orange
graph) and the internal Vensim load model (blue graph).
Discussion
As shown by the figures, the two modeling approaches generates different results on the
modeled utility’s net income, power utilization rate and electricity usage. As can be seen in
Figure 4 and 7 the MATLAB generated load profiles results in larger peak load for a specific
electricity usage. With the increase in peak load the utility has to match with increasing its
generation capacity in order to avoid power utilization rate deterioration. The external
MATLAB load model therefore requires a larger installed capacity to supply the same
electricity consumption, which could have negative impacts on net income. But since total
electricity usage is higher in the MATLAB load model case the utility also receives more
income.
As can be seen in Figure 7, the MATLAB load model results in higher total electricity usage
than the Vensim load model. This might seem unintuitive since the MATLAB load model
requires a higher peak load compared to electricity usage, which therefore requires a larger
generation capacity (increasing expenses that the utility could have otherwise used for
connecting users thereby obtaining a larger total electricity usage) for a given electricity
consumption. However, use of the bottom-up appliance diffusion approach increases
electricity usage faster (as a function of income) than the Vensim approach. In Vensim
increase in electricity usage is proportional to increase in income, e.g. if income increase
with 10% electricity usage increase with 10%. One issue with the Vensim approach is that it
does not take into account acquisition of new appliances, thereby excluding a feedback loop.
The result of the two approaches is a difference in electricity usage sensitivity to income,
which as income grows becomes apparent and can be seen in Figure 7. The increase in
10
electricity usage also increase income for the utility allowing it to make the necessary
investments in new generation.
When comparing the two models it is important to note that the MATLAB load model is also
a simplified load model and excludes certain processes that are common in more advanced
load models. One such processes is coincidence. Coincidence, or coincidence factor, is an
indicator of the probability that loads will be turned on simultaneously. With a coincidence
factor of 0.5, a 50 kW capacity (transmission and generation) is enough to supply an installed
load of 100 kW. As described in the method sections, the standard load profiles used in the
modeling procedure are based on interviews and measurements. Since the measurements
where done at household level the model takes coincidence between appliance loads into
account. However, as the total load profile is aggregated from the individual users, the
model lacks an integration of coincidence between users. If a full coincidence would be
integrated it could allow for more users given the same installed capacity.
Another simplification in the integration between the MATLAB load model and the system
dynamics model is through the definition and implementation of power utilization rate. Due
to the definition of power utilization rate, it only reflects one point in time. Even though this
implementation includes issues with potential power outages it assumes all users are equally
affected. This can be especially troublesome if the peak load that deteriorate power
utilization rate takes place during the evening, when the load is almost exclusively based on
household’s electricity consumption. With a decreased power utilization rate, it will impact
the electricity usage of both households and SMEs, even though SMEs consumption is
mostly during the day and are therefore not affected by the issues to the same extent. These
issues could potentially be solved by implementing different power utilization rate indicators
for each user group, e.g. one for households and one for SMEs.
Conclusions and future work
By incorporating a comparably detailed load model into an existing system dynamics model,
this work has compared the economic impact on a minigrid utility between two different
load modeling methods. The results show that the more detailed load model requires the
utility to invest in generation capacity at more frequent intervals in order to keep power
utilization rate from deteriorating. The results also show that investments in new generation
capacity are more frequent, but that each investment is smaller, when the MATLAB load
model is used. The bottom-up approach to modeling load appliances (MATLAB load model)
also shows a faster growth in electricity consumption
As presented in this work, the improvement of load modeling detail impacts the modelled
minigrid utility’s balance sheet and its investments in new capacity. However, when
matching generation with load it is also important to consider the specific characteristics of
different energy sources, in particular intermittent energy sources. While this has been
excluded in this work, future work could include integration of capacity expansion choices
including renewable energy source characteristics.
11
References
Ahlborg, H., & Hammar, L. (2014). Drivers and barriers to rural electrification in Tanzania and
Mozambique — Grid-extension, off-grid, and renewable energy technologies.
Renewable Energy, 61(0), 117-124. doi:
http://dx.doi.org/10.1016/j.renene.2012.09.057
Azimoh, C. L., Klintenberg, P., Wallin, F., Karlsson, B., & Mbohwa, C. (2016). Electricity for
development: Mini-grid solution for rural electrification in South Africa. Energy
Conversion and Management, 110, 268-277. doi:
http://dx.doi.org/10.1016/j.enconman.2015.12.015
Barnes, D. (2007). The challenge of rural electrification: strategies for developing countries.
Washington, DC: Routledge.
Barnes, D., & Foley, G. (2004). Rural Electrification in the Developing World: A Summary of
Lesson from Successful Programs. [Report]. ESMAP, World Bank.
Bhattacharyya, S. (2013). Rural Electrification Through Decentralised Off-grid Systems in
Developing Countries: Springer.
Bucher, C., & Andersson, G. (2012). Generation of domestic load profiles-an adaptive top-
down approach. Proceedings of PMAPS 2012, 10-14.
Carpinteiro, O. A. S., Leme, R. C., de Souza, A. C. Z., Pinheiro, C. A. M., & Moreira, E. M.
(2007). Long-term load forecasting via a hierarchical neural model with time
integrators. Electric Power Systems Research, 77(3-4), 371-378. doi:
http://dx.doi.org/10.1016/j.epsr.2006.03.014
Davidsen, P. (1992). The structure-behavior diagram: understanding the relationship
between structure and behavior in complex dynamic systems. Paper presented at the
International Conference of the System Dynamics Society, Utrecht, Netherlands.
Diaz, P., Arias, C. A., Pefia, R., & Sandoval, D. (2010). FAR from the grid: A rural electrification
field study. Renewable Energy, - 35(- 12), - 2834.
Goldemberg, J., Johansson, T. B., Reddy, A. K. N., & Williams, R. H. (1985). Basic Needs and
Much More with One Kilowatt per Capita. Ambio, 14(4/5), 190-200. doi:
10.2307/4313148
Hartvigsson, E. (2016). A System Dynamics Analysis of Cost-Recovery: A Study of Rural
Minigrid Utilities in Tanzania. Licentiate, Chalmers.
Hartvigsson, E., Ehnberg, J., Ahlgren, E., & Molander, S. (2015). Assessment of load profiles in
minigrids: A case in Tanzania. Paper presented at the Power Engineering Conference
(UPEC), 2015 50th International Universities.
Hjorth, P., & Bagheri, A. (2006). Navigating towards sustainable development: A system
dynamics approach. Futures, 38(1), 74-92.
IEA. (2013). World Energy Outlook: International Energy Agency.
Kirubi, C., Jacobson, A., Kammen, D. M., & Mills, A. (2009). Community-Based Electric Micro-
Grids Can Contribute to Rural Development: Evidence from Kenya. World
12
Development, 37(7), 1208-1221. doi:
http://dx.doi.org/10.1016/j.worlddev.2008.11.005
Levin, T., & Thomas, V. M. (2014). Utility-maximizing financial contracts for distributed rural
electrification. Energy, 69(0), 613-621. doi:
http://dx.doi.org/10.1016/j.energy.2014.03.057
Madubansi, M., & Shackleton, C. M. (2006). Changing energy profiles and consumption
patterns following electrification in five rural villages, South Africa. Energy Policy,
34(18), 4081-4092. doi: 10.1016/j.enpol.2005.10.011
Marzouk, M., & Azab, S. (2014). Environmental and economic impact assessment of
construction and demolition waste disposal using system dynamics. Resources,
Conservation and Recycling, 82, 41-49. doi:
http://dx.doi.org/10.1016/j.resconrec.2013.10.015
McNeil, M. A. (2008). Global potential of energy efficiency standards and labeling programs.
Lawrence Berkeley National Laboratory.
Niu, D., Li, J., Li, J., & Liu, D. (2009). Middle-long power load forecasting based on particle
swarm optimization. Computers & Mathematics with Applications, 57(11-12), 1883-
1889. doi: http://dx.doi.org/10.1016/j.camwa.2008.10.044
Qudrat-Ullah, H., & Davidsen, P. |. (2001). Understanding the dynamics of electricity supply,
resources and pollution: Pakistan's case. Energy, 26(6), 595-606. doi:
http://dx.doi.org/10.1016/S0360-5442(01)00019-6
Sarangi, G. K., Pugazenthi, D., Mishra, A., Palit, D., Kishore, V. V. N., & Bhattacharyya, S. C.
(2014). Poverty Amidst Plenty: Renewable Energy-Based Mini-Grid Electrification in
Nepal. In S. C. Bhattacharyya & D. Palit (Eds.), Mini-Grids for Rural Electrification of
Developing Countries: Analysis and Case Studies from South Asia (pp. 343-371).
Schnitzer, D., Lounsbury S., D., Carvallo, J. P., Deshmukh, R., Apt, J., & Kammen M., D. (2014).
Microgrids for Rural Electrification: A critical review of best practices based on seven
case studies: United Nations Foundation.
Stokes, M. (2005). Removing barriers to embedded generation: a fine-grained load model to
support low voltage network performance analysis. PhD, De Montfort University.
Tenenbaum, B., Greacen, C., Siyambalapitya, T., & Knuckles, J. (2014). From the Bottom Up -
How Small Power Producers and Mini-grids Can Deliver Electrification and Renewable
Energy in Africa. [Book]. The World Bank.
Teufel, F., Miller, M., Genoese, M., & Fichtner, W. (2013). Review of System Dynamics
models for electricity market simulations.
Urpelainen, J. (2014). Grid and off-grid electrification: An integrated model with applications
to India. Energy for Sustainable Development, 19, 66- 71.
Widén, J., Lundh, M., Vassileva, |., Dahlquist, E., Ellegard, K., & Wackelgard, E. (2009).
Constructing load profiles for household electricity and hot water from time-use
data—Modelling approach and validation. Energy and Buildings, 41(7), 753-768. doi:
http://dx.doi.org/10.1016/j.enbuild.2009.02.013
13
of the simulation experience. Having other students with whom to discuss their ideas
fostered an interactive, engaging experience to test ideas in a non-threatening, low-stakes
environment. The experiences of the students who have used the simulation thus far and
the content of this interview with just a few students is just a small sampling of the
potential for learning.
Certainly, additional data will be useful in determining whether these reflections are truly
representative of teachers and students in general. However, going back to the beginnings
of this work in K-12, we know that students tend to find these types of experiences
interesting and empowering. These middle school science teachers who used similar
simulations starting more than 25 years ago state, “Students shift from being passive
consumers of information that need only be remembered until the next exam, to being
active participants in the acquisition and utilization of knowledge. In this environment
they become intrinsically motivated to extend the information they have learned” (Draper
and Swanson, 1990).
Conclusion
After adapting the model, building a new interface, testing within classrooms in both
middle and high school and developing a teacher guide with handouts, this new resource is
available for widespread use. Educators can access the simulation and download the free
guide from the Creative Learning Exchange (CLE) website:
http://www.clexchange.org/curriculum/simulations/prison_simulation.asp.
Because of the highly engaging, interactive, and self-directed nature of simulations such as
this one, a challenge is to continue to develop other useful models with accompanying
experiences that align with K-12 curricular goals. This challenge would be best met
through identifying future partnerships with individuals or groups who could fund such
efforts, experienced Systems Dynamics modelers, K-12 teachers and K-12 implementers of
systems thinking/dynamic modeling. The team assembled for this project brought
together that exact combination of experience.
Page 14 of 21
Barry also created a model (Richmond, 1977) for the Milgram Experiments (Milgram,
1963) that has strong connections to understanding historical events and themes studied
in K-12 and would also create an engaging context for students to explore. As one student
states in the interview, “It’s so much more powerful to be able to experience and play
around with these factors yourself and predict what will happen, because it really gets the
creative juices flowing, rather than just reading about it. These can definitely apply to real
life - you see corruption and fear all the time in history.” Now is the time to consider, how
can others share their expertise as Barry Richmond chose to do over many years within K-
12? How can others like Barry, learn from the experiences and expertise of K-12
educators? How can we together continue to create similar valuable learning experiences
for this generation of youth?
Page 15 of 21
Supporting Materials
Behind Closed Gates : CC2014_Behind Closed Gates.pdf
Behind Closed Gates Simulation Model: BCG-modelonly.stmx
The entire Behind Closed Gates curriculum is available online:
http://www.clexchange.org/curriculum/simulations/prison_simulation.asp
Page 16 of 21
Appendix A: Revised Model Equations
{VERSION 10.0.4 }
{ INITIALIZATION EQUATIONS }
:c Initial_guard_distrust = 0.1
: Ss Guard_distrust = Initial_guard_distrust
:cInitial_guard_repression = 0.1
: sGuard_repression = Initial_guard_repression
: c Initial_prisoner_fear = 0.1
: S Prisoner_fear = Initial_prisoner_fear
: c Initial_prisoner_resistance = 0.1
: S Prisoner_resistance = Initial_prisoner_resistance
: c Initial_prisoner_solidarity = 0.1
: s Prisoner_solidarity = Initial_prisoner_solidarity
: cNormal_guard_perception_time = 0.25
: c Impact_of_resistance_on_perception_time = GRAPH(Prisoner_resistance / ( Guard_distrust + 0.001 ))
(0.9, 120), (0.92, 120), (0.94, 110.933333333), (0.96, 88), (0.98, 40), (1, 855555555556), (1.02, 1), (1.04, 1),
(1.06, 1), (1.08, 1), (1.1, 1)
: cGuard__perception_time = Normal_guard__perception_time * Impact_of_resistance_on_perception_time
: c Tendency_to_distrust = 1
: c Impact_of_unity_on_distrust = GRAPH(Prisoner_solidarity)
(0, 0.75), (0.1, 0.761904761905), (0.2, 0.792857142857), (0.3, 0.852380952381), (0.4, 0.969047619048),
(0.5, 1.20476190476), (0.6, 1.36428571429), (0.7, 1.42857142857), (0.8, 1.46904761905), (0.9,
149047619048), (1, 1.5)
: b Change_in_guard_distrust = ( (( Prisoner_resistance * Tendency_to_distrust ) - Guard_distrust) *
Impact_of_unity_on_distrust )/ Guard_perception_time
: c Impact_of_unity_on_fear = GRAPH(Prisoner_solidarity)
(0, 1), (0.1, 0.979), (0.2, 0.958), (0.3, 0.932), (0.4, 0.903), (0.5, 0.866), (0.6, 0.831), (0.7, 0.785), (0.8, 0.76),
(0.9, 0.75), (1, 0.75)
: c Indicated_fear_from_repression = GRAPH(Guard_repression)
(0, 0.00952380952381), (0.1, 0.02), (0.2, 0.045), (0.3, 0.089), (0.4, 0.179), (0.5, 0.301), (0.6, 0.459), (0.7,
0.626), (0.8, 0.749), (0.9, 0.863), (1, 0.971428571429)
: c Indicated__prisoner_fear = Indicated_fear_from_repression * Impact_of_unity__on_fear
: c Tendency_to_fear = 1
: ¢ Time_to_change_fear = 0.5
: b Change_in__prisoner_fear = ( ( Indicated__prisoner_fear * Tendency_to_fear ) - Prisoner_fear) /
Time_to_change_fear
: c Indicated_level_of_repression_by_guards = GRAPH(Guard_distrust)
(0, 0), (0.1, 0.053), (0.2, 0.122), (0.3, 0.256), (0.4, 0.439), (0.5, 0.65), (0.6, 0.772), (0.7, 0.878), (0.8, 0.939),
(0.9, 0.976), (1, 1)
: c Time_to_change_repressiveness = 0.5
: c Willingness_to_repress = 1
: b Change_in_repression = ( ( Indicated_level_of_repression_by_guards * Willingness_to_repress ) -
Guard_repression ) / Time_to_change__repressiveness
: c Fraction_decrease_from_fear = GRAPH(Prisoner_fear)
(0, 0), (0.1, 0.02), (0.2, 0.065), (0.3, 0.23), (0.4, 0.445), (0.5, 1), (0.6, 1.3275), (0.7, 1.5), (0.8, 1.5), (0.9, 1.5), (1,
15)
: c Fraction_decrease_from_repression = GRAPH(Guard_repression)
(0, 0), (0.1, 0.01), (0.2, 0.01), (0.3, 0.02), (0.4, 0.03), (0.5, 0.08), (0.6, 0.24), (0.7, 2.835), (0.8, 3), (0.9, 3), (1,3)
: ¢ Fraction_decrease_in_resistance = Fraction_decrease_from_fear + Fraction_decrease_from_repression
: f Decrease_in_resistance = Prisoner_resistance * Fraction_decrease_in_resistance
: c Impact_of_fear_on_unity_decrease = GRAPH(Prisoner_fear)
(0, 1), (0.1, 1.472), (0.2, 1.825), (0.3, 2.65), (0.4, 3.947), (0.5, 5.833), (0.6, 10.667), (0.7, 19.321), (0.8, 26.133),
(0.9, 29.263), (1, 30)
:¢ Time_for_unity__to_dissolve = 3
Page 17 of 21
: fDecrease_in_unity = ( Prisoner_solidarity / Time_for_unity__to_dissolve ) *
Impact_of_fear_on_unity_decrease
: c Impact_of_fear_on_resistance = GRAPH (Prisoner_fear)
(0, 1), (0.1, 1), (0.2, 0.987301587302), (0.3, 0.946031746032), (0.4, 0.88253968254), (0.5, 0.8), (0.6, 0.65),
(0.7, 0.445), (0.8, 0.275), (0.9, 0.11746031746), (1, 0)
: c Impact_of_repression_on_resistance = GRAPH(Guard_repression)
(0, 0), (0.1, 0.47619047619), (0.2, 0.736507936508), (0.3, 0.869841269841), (0.4, 0.92), (0.5, 0.97), (0.6,
0.985), (0.7, 1), (0.8, 1), (0.9, 1), (1, 1)
: c Indicated__resistance = Impact_of_repression__on_resistance * Impact_of_fear_on_resistance
:c Time_to_change_resistance = 1
: c Willingness_to_resist = 1
: f Increase__in_resistance = ( ( Indicated_resistance * Willingness_to_resist ) - Prisoner_resistance ) /
Time_to_change_resistance
: c Indicated_unity_from_repression = GRAPH(Guard_repression)
(0, 0.87), (0.1, 0.813), (0.2, 0.72), (0.3, 0.573), (0.4, 0.374), (0.5, 0.24), (0.6, 0.15), (0.7, 0.098), (0.8, 0.053),
(0.9, 0.012), (1, 0)
: c Tendency_to_stand_together = 1
: ¢ Time_to_change_unity = 1
: fIncrease__in_unity = ( ( Indicated_unity__from_repression * Tendency_to_stand_together ) -
Prisoner_solidarity ) / Time_to_change_unity
: c Fear_prediction = 0
: c Guard_Distrust_Prediction = GRAPH(TIME)
(0, -1), (1.5, -1), (3, -1), (4.5, -1), (6, -1), (7.5, -1), (9, -1), (10.5, -1), (12, -1), (13.5, -1), (15, -1)
: cGuard_Repression_Prediction = GRAPH(TIME)
(0, -1), (1.5, -1), (3, -1), (4.5, -1), (6, -1), (7.5, -1), (9, -1), (10.5, -1), (12, -1), (13.5, -1), (15, -1)
: c Prisoner_Fear_Prediction = GRAPH(TIME)
(0, -1), (1.25, -1), (2.5, -1), (3.75, -1), (5, -1), (6.25, -1), (7.5, -1), (8.75, -1), (10, -1), (11.25, -1), (12.5, -1),
(13.75, -1), (15, -1)
:c Prisoner_Resistance_Prediction = GRAPH(TIME)
(0, -1), (1.5, -1), (3, -1), (4.5, -1), (6, -1), (7.5, -1), (9, -1), (10.5, -1), (12, -1), (13.5, -1), (15, -1)
: c Prisoner_Solidarity_Prediction = GRAPH(TIME)
(0, -1), (1.5, -1), (3, -1), (4.5, -1), (6, -1), (7.5, -1), (9, -1), (10.5, -1), (12, -1), (13.5, -1), (15, -1)
{ RUNTIME EQUATIONS }
:s Guard_distrust(t) = Guard_distrust(t - dt) + (Change_in_guard_distrust) * dt
: s Guard_repression(t) = Guard_repression(t - dt) + (Change_in_repression) * dt
: s Prisoner_fear(t) = Prisoner_fear(t - dt) + (Change_in_prisoner_fear) * dt
: s Prisoner_resistance(t) = Prisoner_resistance(t - dt) + (Increase_in_resistance - Decrease_in_resistance) *
dt
: s Prisoner_solidarity(t) = Prisoner_solidarity(t - dt) + (Increase_in_unity - Decrease_in_unity) * dt
: c Impact_of_resistance_on_perception_time = GRAPH(Prisoner_resistance / ( Guard_distrust + 0.001 ))
(0.9, 120), (0.92, 120), (0.94, 110.933333333), (0.96, 88), (0.98, 40), (1, 855555555556), (1.02, 1), (1.04, 1),
(1.06, 1), (1.08, 1), (1.1, 1)
:c Guard_perception_time = Normal_guard_perception_time * Impact_of_resistance_on_perception_time
: c Impact_of_unity__on_distrust = GRAPH(Prisoner__solidarity)
(0, 0.75), (0.1, 0.761904761905), (0.2, 0.792857142857), (0.3, 0.852380952381), (0.4, 0.969047619048),
(0.5, 1.20476190476), (0.6, 1.36428571429), (0.7, 1.42857142857), (0.8, 1.46904761905), (0.9,
1.49047619048), (1, 1.5)
: b Change_in_guard_distrust = ( (( Prisoner_resistance * Tendency_to_distrust ) - Guard_distrust) *
Impact_of_unity_on_distrust )/ Guard_perception_time
: c Impact_of_unity__on_fear = GRAPH(Prisoner_solidarity)
(0, 1), (0.1, 0.979), (0.2, 0.958), (0.3, 0.932), (0.4, 0.903), (0.5, 0.866), (0.6, 0.831), (0.7, 0.785), (0.8, 0.76),
(0.9, 0.75), (1, 0.75)
:c Indicated_fear_from_repression = GRAPH(Guard_repression)
Page 18 of 21
(0, 0.00952380952381), (0.1, 0.02), (0.2, 0.045), (0.3, 0.089), (0.4, 0.179), (0.5, 0.301), (0.6, 0.459), (0.7,
0.626), (0.8, 0.749), (0.9, 0.863), (1, 0.971428571429)
: c Indicated__prisoner_fear = Indicated_fear_from_repression * Impact_of_unity__on_fear
: b Change_in__prisoner_fear = ( ( Indicated_prisoner_fear * Tendency_to_fear ) - Prisoner_fear) /
Time_to_change_fear
: c Indicated_level_of_repression_by_guards = GRAPH(Guard_distrust)
(0, 0), (0.1, 0.053), (0.2, 0.122), (0.3, 0.256), (0.4, 0.439), (0.5, 0.65), (0.6, 0.772), (0.7, 0.878), (0.8, 0.939),
(0.9, 0.976), (1, 1)
: b Change_in_repression = ( ( Indicated_level_of_repression_by_guards * Willingness_to_repress ) -
Guard_repression ) / Time_to_change_repressiveness
: c Fraction_decrease_from_fear = GRAPH(Prisoner_fear)
(0, 0), (0.1, 0.02), (0.2, 0.065), (0.3, 0.23), (0.4, 0.445), (0.5, 1), (0.6, 1.3275), (0.7, 1.5), (0.8, 1.5), (0.9, 1.5), (1,
15)
: c Fraction_decrease_from_repression = GRAPH(Guard_repression)
(0, 0), (0.1, 0.01), (0.2, 0.01), (0.3, 0.02), (0.4, 0.03), (0.5, 0.08), (0.6, 0.24), (0.7, 2.835), (0.8, 3), (0.9, 3), (1, 3)
: c Fraction_decrease_in_resistance = Fraction_decrease_from_fear + Fraction_decrease_from_repression
: fDecrease_in_resistance = Prisoner_resistance * Fraction_decrease_in_resistance
: c Impact_of_fear_on_unity_decrease = GRAPH(Prisoner_fear)
(0, 1), (0.1, 1.472), (0.2, 1.825), (0.3, 2.65), (0.4, 3.947), (0.5, 5.833), (0.6, 10.667), (0.7, 19.321), (0.8, 26.133),
(0.9, 29.263), (1, 30)
: fDecrease_in_unity = ( Prisoner_solidarity / Time_for_unity__to_dissolve ) *
Impact_of_fear__on_unity_decrease
: c Impact_of_fear_on_resistance = GRAPH(Prisoner_fear)
(0, 1), (0.1, 1), (0.2, 0.987301587302), (0.3, 0.946031746032), (0.4, 0.88253968254), (0.5, 0.8), (0.6, 0.65),
(0.7, 0.445), (0.8, 0.275), (0.9, 0.11746031746), (1, 0)
: c Impact_of_repression__on_resistance = GRAPH (Guard_repression)
(0, 0), (0.1, 0.47619047619), (0.2, 0.736507936508), (0.3, 0.869841269841), (0.4, 0.92), (0.5, 0.97), (0.6,
0.985), (0.7, 1), (0.8, 1), (0.9, 1), (1, 1)
: c Indicated_resistance = Impact_of_repression__on_resistance * Impact_of_fear_on_resistance
: fIncrease_in_resistance = ( ( Indicated__resistance * Willingness_to_resist ) - Prisoner_resistance ) /
Time_to_change_resistance
: c Indicated_unity_from_repression = GRAPH(Guard_repression)
(0, 0.87), (0.1, 0.813), (0.2, 0.72), (0.3, 0.573), (0.4, 0.374), (0.5, 0.24), (0.6, 0.15), (0.7, 0.098), (0.8, 0.053),
(0.9, 0.012), (1, 0)
: fIncrease__in_unity = ( ( Indicated_unity_from_repression * Tendency_to_stand_together ) -
Prisoner_solidarity ) / Time_to_change_unity
: cGuard_Distrust_Prediction = GRAPH(TIME)
(0, -1), (1.5, -1), (3, -1), (4.5, -1), (6, -1), (7.5, -1), (9, -1), (10.5, -1), (12, -1), (13.5, -1), (15, -1)
: c Guard_Repression_Prediction = GRAPH(TIME)
(0, -1), (1.5, -1), (3, -1), (4.5, -1), (6, -1), (7.5, -1), (9, -1), (10.5, -1), (12, -1), (13.5, -1), (15, -1)
: c Prisoner_Fear_Prediction = GRAPH(TIME)
(0, -1), (1.25, -1), (2.5, -1), (3.75, -1), (5, -1), (6.25, -1), (7.5, -1), (8.75, -1), (10, -1), (11.25, -1), (12.5, -1),
(13.75, -1), (15, -1)
: c Prisoner_Resistance_Prediction = GRAPH(TIME)
(0, -1), (1.5, -1), (3, -1), (4.5, -1), (6, -1), (7.5, -1), (9, -1), (10.5, -1), (12, -1), (13.5, -1), (15, -1)
:c Prisoner_Solidarity_Prediction = GRAPH(TIME)
(0, -1), (1.5, -1), (3, -1), (4.5, -1), (6, -1), (7.5, -1), (9, -1), (10.5, -1), (12, -1), (13.5, -1), (15, -1)
{ TIME SPECS }
STARTTIME=0
STOPTIME=14
DT=0.03125
INTEGRATION=RK4
RUNMODE=NORMAL
PAUSEINTERVAL=INF
Page 19 of 21
Appendix B: Example Student Test Run and Response
Note: Predictions on the graphs are in red. Student responses are in italics.
Run 3:
1. Variable Changed: Fear
2. Setting: 0.0
3. Prediction: How do you think the five variables will change over time compared to the
baseline run? Draw your predictions on the prediction screen or describe them here.
Because the fear is down they won't be afraid to resist, so resistance and solidarity will go up,
but because they do distrust, repression will go up.
Distrust Repressio Fear Solidari Resistance
Ale
4. After running, describe what actually happened in comparison to the base run. [insert a copy
of the graphs, either on paper or electronically.]
Distrust and repression were close to what I expected, but the solidarity and resistance shoots
up, then falls.
In the control room, you set the five sliders
and run the simulation over the 14 days.
Make sure to make predictions on the
predictions screen before you hit the run
The graphs will show what actually happens
in blue and your predictions in red
5. Why did this happen?
I get why distrust and repression go up and slowly go down, but I keep getting the resistance and
solidarity wrong. I probably am not thinking about that when they go up and the distrust and
repression go up, they are going to fall.
Page 20 of 21
References
Andersen, J., LaVigne, A, Potash, J, and Stuntz, L.(2014). Behind Closed Gates Simulation and
Lesson Guide. Creative Learning Exchange.
Collins, Suzanne (1996). The Hunger Games. Scholastic.
Doyle, James with Khalid Saeed and Jeanine Skorinko (2009). Personal versus Situational
Dynamics: Implications of Barry Richmond's Models of Classic Experiments in Social
Psychology. System Dynamics Society.
Draper, F. and Swanson, M. (1990). Learner-directed Systems Education: A Successful
Example. System Dynamics Review, 6(2), 209-213
Forrester, Jay (2009). Learning through System Dynamics as Preparation for the 21st
Century (D-4895), Creative Learning Exchange
(http://www.clexchange.org/ftp/documents/whyk12sd/Y_2009-
02LearningThroughSD.pdf)
Golding, William (1963). Lord of the Flies. Faber & Faber.
Milgram, S. (1963). Behavioral Study of Obedience. Journal of Abnormal and Social
Psychology, 67(4):371-378.
Orwell, George (1946). Animal Farm. Harcourt, Brace and Company. New York.
Richmond, B. (1977). Generalization with Individual Uniqueness: Modeling the Milgram
Experiments. System Dynamics Group.
Richmond, B. (1993). Systems Thinking: Critical Thinking Skills for the 1990s and Beyond.
System Dynamics Review 9(2), 113-133.
Zimbardo, P. G., Haney, C., Banks, C., & Jaffe, D. (1973). The Mind is a Formidable Jailer: A
Pirandellian Prison. The New York Times Magazine, pp. 38 ff.
Zimbardo, P.G. (2005). The Psychology of Power and Evil: All Power to the Person? To the
Situation? To the System?. Manuscript, Psychology Department, Stanford University.
Page 21 of 21
