---
title: "System Dynamics"
author: "Ricardo R. Palma"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


System Dynamics
CPT Rick Mard
August 30, 2017

## 1 Introduction
This tutorial is intended to introduce system dynamics modeling concepts and techniques using the R programming language. The concepts and techniques presented herewithin are simple examples designed for the reader to follow along using R Studio. Therefore, both R and R Studio are suggested in order to fully benefit from the proceeding tutorial.

### 1.1 Packages
The following packages are used throughout the tutorial. The deSolve package is primarily used to calculate the differential equations and numerical integration. Other required packages include ggplot2, gridextra, plyr, magrittr, and grid. These packages support visualizations, data frame manipulation, and pipe operators.

```{r}
library(deSolve) #supports numerical integration using a range of numerical methods
library(ggplot2) #supports visualization of layered graphics
library(gridExtra) #supports visualization of multiple plots & graphs
library(plyr) #supports data frame merging
library(magrittr) #supports pipe function
library(grid) #supports additional graphing capability
```



## 2 Background de System Dynamics 

La dinámica de sistemas es una metodología aplicable al modelado de factores de un sistema complicado de factores interdependientes como el que se presenta en este tutorial. Jay Forrester descubrió la dinámica de sistemas después de su trabajo con la gerencia de General Electric y con sistemas de radar para la Marina de los EE. UU. Combinó las complejidades de los sistemas sociales con los sistemas físicos en la recientemente creada Escuela de Administración del Instituto Tecnológico de Massachusetts, o MIT (Forrester 1968). Forrester describe los sistemas como

una agrupación de partes que operan juntas para un propósito común (System Dynamics Society 2017).

El Departamento de Defensa opera de manera jerárquica con información y conocimiento controlados por clasificación y necesidad. En relación con esto, Forrester informa que “sin una estructura organizativa, el conocimiento es una mera colección de observaciones, prácticas e incidentes conflictivos” (System Dynamics Society 2017). Para llegar a soluciones para los esfuerzos de planificación estratégica y operativa, se necesitan sistemas que definan las complejidades y relaciones dentro del entorno operativo.
### 2.1 Models
A system dynamics model contains causal loops that identify the feedback of interactions between elements of a system. These causal loops are usually converted to ‘stocks’ and ‘flows’. Stocks represent elements with capacities that change over time. Flows define the rate at which a stock can change.

Stock and flow diagrams display the interactions between elements in the system so information can be traced throughout. Linear and/or continuous equations in time are applied to the system flows before initializing the system with appropriate stock levels based on expert input, previous research, or analyses (Duggan 2016).

### 2.2 Stocks & Flows

Un stock es la base de cualquier sistema (Meadows 2008), y los stocks y los flujos son los componentes básicos de los modelos de dinámica de sistemas. Caracterizan el estado del sistema bajo estudio, además de proporcionar la información en la que se basan las decisiones y acciones (Sterman 2000). Las acciones solo pueden cambiar a través de sus flujos, que son las cantidades agregadas (entradas) o restadas (salidas) de una acción a lo largo del tiempo (Duggan 2016).

![Stock and flow model of customers](fig1.png)

    
    Stock and flow model of customers

### 2.3 Feedback

Un stock es la base de cualquier sistema (Meadows 2008), y los stocks y los flujos son los componentes básicos de los modelos de dinámica de sistemas. Caracterizan el estado del sistema bajo estudio, además de proporcionar la información en la que se basan las decisiones y acciones (Sterman 2000). Las acciones solo pueden cambiar a través de sus flujos, que son las cantidades agregadas (entradas) o restadas (salidas) de una acción a lo largo del tiempo (Duggan 2016).
## 3 deSolve Package

El paquete deSolve de R resuelve problemas de valores iniciales escritos como ecuaciones diferenciales ordinarias (ODE), ecuaciones algebraicas diferenciales (DAE) y ecuaciones diferenciales parciales (PDE) Soetaert et al. (2010). Para los modelos de dinámica de sistemas, se utiliza el solucionador ODE en deSolve. El requisito clave es que los modeladores de dinámica de sistemas implementen las ecuaciones del modelo en una función, y esta función es llamada por deSolve (Duggan 2016).
Primero, inicializamos las constantes de tiempo de simulación y creamos un vector de tiempo de simulación usando seq():

```{r}
# Set time period and step
START <- 2020
FINISH <- 2050
STEP <- 0.25
simtime <- seq(START, FINISH, by = STEP)

head(simtime)
```



```{r}
tail(simtime)
```



Los vectores del modelo incluyen acciones y su valor inicial y eje, que son los parámetros exógenos.


```{r}
# Set stock capacity, growth and decline rates
stocks <- c(sCustomers = 10000)
auxs <- c(aGrowthFraction = 0.08, aDeclineFraction = 0.03)
```


Se crea una función para ejecutar las ecuaciones para el tiempo de simulación, el vector de stock actual y los parámetros del modelo. Estos se almacenan con as.list() para acceder. Se trata de flujos de contrataciones y pérdidas junto con el flujo neto, un derivado. El marco de datos 'o' se define para permitir que la función ode del paquete deSolve resuelva las ecuaciones anteriores. Los resultados resultantes incluyen estadísticas resumidas para visualizaciones.

```{r}
# Create the function and returns a list
model <- function(time, stocks, auxs){
  with(as.list(c(stocks, auxs)), {
    fRecruits <- sCustomers * aGrowthFraction
    fLosses <- sCustomers * aDeclineFraction
    dC_dt <- fRecruits - fLosses
    return (list(c(dC_dt),
                 Recruits = fRecruits, Losses = fLosses,
                 GF = aGrowthFraction, DF = aDeclineFraction))
  }) 
}

# Create the data frame using the `ode` function
o <- data.frame(ode(y = stocks, times = simtime, func = model,
                    parms = auxs, method = "euler"))

head(o)
```


```{r}
tail(o)
summary(o[,-c(1, 5, 6)])
```



En el siguiente gráfico, el número de clientes es el resultado de la ecuación diferencial de contrataciones y pérdidas definida en el modelo de función anterior.


### Plot the output

```{r}
o %>% 
  ggplot() +
  geom_line(aes(time, o$sCustomers), colour = "blue") +
  geom_point(aes(time, o$sCustomers), colour = "blue") +
  scale_y_continuous() +
  ylab("Customers") +
  xlab("Year") +
  ggtitle("Customers Over Time", subtitle = "2020 - 2050")

```



## 4 Limits to Growth

Es necesario introducir un método de formulación antes de presentar los límites del modelo de crecimiento. Los modelos de formulación permiten a los desarrolladores crear ecuaciones lo suficientemente sólidas como para modelar el efecto de una variable sobre otra. Son útiles para una amplia variedad de modelos de dinámica de sistemas, especialmente aquellos con acciones influyentes en diferentes sistemas.

Comenzamos modelando relaciones causales utilizando efectos. Esto se define como

$$Y=Y∗∗Effect(X1onY)∗...∗Effect(XnonY)$$


donde, basado en el supuesto:

* $Y$ = dependent variable of a causal relationship and is a function of n independent variables (X1,X2,...,Xn)

* $Y∗$  = reference value, which is the normal value the variable Y takes on. This reference value is multiplied by a sequence of effect functions that are calculated based on the normalized ratio of input term $(Xi/X∗i)$
 where $X∗i$  is the reference input value, and $Xi$
is the actual input value.

La función de efecto tiene la relación normalizada (X/X∗)
  en su eje x, y siempre contiene el punto (1, 1). Este punto (1, 1) es importante por la siguiente razón: si x
  es igual a su valor de referencia X∗
, entonces la función de efecto será 1, y por lo tanto y
  entonces será igual a su valor de referencia Y∗
.

### 4.1 Example




Usando un valor de X asignado (borde),
$$GrowthRate=RefGrowthRate∗EffectofAvailabilityonGrowthRate$$

$$EffectofAvailabilityonGrowthRate=f(AvailabilityRefAvailability)$$

* RefGrowthRate=0.10
* RefAvailability=1.0


Relationship between availability and effect on growth rate
Relationship between availability and effect on growth rate


La relación lineal entre disponibilidad y tasa de crecimiento en los dos casos extremos anteriores es:

* Availability = 1, Effect = 1, so growth rate is maximized
* Availability = 0, Effect = 0, since there are no resources, therefore growth rate is zero

which keeps the model simple and follows the slope-intercept equation of $y=mx+c$
. The effect equation can be modified from above.



$$EffectofAvailabilityonGrowthRate=AvailabilityRefAvailability$$


La tasa de crecimiento depende del valor de la disponibilidad a través del efecto variable, como se ve en la siguiente tabla.

Ref availability |	Availability|	Effect of availability on growth rate	|Ref growth rate	|Growth rate
-----------------|--------------|---------------------------------------|----------------|
1.0	|1.0	|1.0	|0.10	|0.10
1.0	|0.5	|0.5	|0.10	|0.5
1.0	|0.0	|0.0	|0.10	|0.0


### 4.2 S-Shaped Growth

Del modelo anterior, se agrega un factor limitante como un circuito de equilibrio para contrarrestar el crecimiento como en el modelo siguiente

### 4.2.1 Define Variables

[A one-stock model of limits to growth](fig2.png)




donde la acción tiene una entrada con un valor inicial de 100

$Stock=\int_{100}^n NetFlow$


y el flujo de entrada 

$$NetFlow=Stock∗GrowthRate$$


con capacidad del sistema relación de disponibilidad

$$Availability=1−\frac{Stock}{Capacity}$$


and a capacity of 10000 (arbitrary) to limit growth when stock equals capacity since availability will be zero

    Capacity=10,000


### 4.2.2 Declare Variables using deSolve

Using the three R libraries, declare stock and auxiliary variables for the simulation including time period and step.

```{r}
library(deSolve)
library(ggplot2)
library(gridExtra)
```


```{r}
START <- 0
FINISH <- 100
STEP <- 0.25
simtime <- seq(START, FINISH, by = STEP)
stocks <- c(sStock = 100)
auxs <- c(aCapacity = 10000, 
          aRef.Availability = 1,
          aRef.GrowthRate = 0.10)
```


Set the time period and step. Define the stocks and auxiliaries.

### 4.2.3 Build a Function

The model’s equations are contained within a function called by deSolve by timestep and passes a vector of stocks with sim values and a vector of auxiliaries. These values are evaluated from the equations (above). The variable dS_dt is an integral. The data frame sModel is created and the first few lines are displayed.

```{r}
# Create the function
model <- function(time, stocks, auxs){
  with(as.list(c(stocks, auxs)),{
    aAvailability <- 1 - sStock / aCapacity
    aEffect <- aAvailability / aRef.Availability
    aGrowth.Rate <- aRef.GrowthRate * aEffect
    fNet.Flow <- sStock * aGrowth.Rate
    dS_dt <- fNet.Flow
    return(list(c(dS_dt), NetFlow = fNet.Flow,
                GrowthRate = aGrowth.Rate, 
                Effect = aEffect,
                Availability = aAvailability))
  })
}
```






     Create the data frame

```{r}
sModel <- data.frame(ode(y = stocks, times = simtime, func = model,
                    parms = auxs, method = "euler"))

head(sModel)
```


```{r}
# Run the following snippet to enable `multiplot`
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
```


### 4.2.4 Visualize the Variables

The data frame sModel is plotted using ggplot2 and a separate multiplot function to display the following relationships. In the Stock versus Time plot, the classic growth mode is exhibited by an exponential growth to time period 46 followed by a logarithmic growth to 100. The Availability plot mirrors the Stock plot since the resource diminishes as stock increases. The Growth Rate is driven by the effect function

$$ EffectofAvailabilityonGrowthRate=f(AvailabilityRefAvailability),$$
$$GrowthRate=RefGrowthRate∗EffectofAvailabilityonGrowthRate$$


and follows Availability until it reaches its capacity. Finally, Net Flow drives the Stock value by rising and falling exponentially with the rate of change.

```{r}
# Plot the outputs using `multiplot`
stockTimePlot <- sModel %>% 
  ggplot() +
  geom_line(aes(time, sStock), color = "purple")

availTimePlot <- sModel %>% 
  ggplot() +
  geom_line(aes(time, Availability), color = "blue")

netflowTimePlot <- sModel %>% 
  ggplot() +
  geom_line(aes(time, NetFlow), color = "red")

growthrateTimePlot <- sModel %>% 
  ggplot() +
  geom_line(aes(time, GrowthRate), color = "green")

multiplot(stockTimePlot, netflowTimePlot, availTimePlot, growthrateTimePlot, cols = 2)


```



The above model is related to the differential equation formulated by Belgian mathematician Verhulst in 1845 (below) to model population increase with growth rate r and limit K:


$$\frac{dP}{dt}=rP∗(1−\frac{P}{K})$$


This model assumes the capacity to be constant, which is problematic because the stock is consumed over time in most realistic models (Duggan 2016).

### 4.2.5 Exercise

What happens when the growth rate is increased to 25%? Show the output in a multiplot.

```{r}
# Change the growth rate to 0.25
auxs <- c(aCapacity = 10000, 
          aRef.Availability = 1,
          aRef.GrowthRate = 0.25)

# Rerun the function
model <- function(time, stocks, auxs){
  with(as.list(c(stocks, auxs)),{
    aAvailability <- 1 - sStock / aCapacity
    aEffect <- aAvailability / aRef.Availability
    aGrowth.Rate <- aRef.GrowthRate * aEffect
    fNet.Flow <- sStock * aGrowth.Rate
    dS_dt <- fNet.Flow
    return(list(c(dS_dt), NetFlow = fNet.Flow,
                GrowthRate = aGrowth.Rate, 
                Effect = aEffect,
                Availability = aAvailability))
  })
}

```

```{r}
# Resave the data frame with the updated auxiliaries and function output
sModel <- data.frame(ode(y = stocks, times = simtime, func = model,
                    parms = auxs, method = "euler"))

# Run `multiplot`
stockTimePlot <- sModel %>% 
  ggplot() +
  geom_line(aes(time, sStock), color = "purple")

availTimePlot <- sModel %>% 
  ggplot() +
  geom_line(aes(time, Availability), color = "blue")

netflowTimePlot <- sModel %>% 
  ggplot() +
  geom_line(aes(time, NetFlow), color = "red")

growthrateTimePlot <- sModel %>% 
  ggplot() +
  geom_line(aes(time, GrowthRate), color = "green")

multiplot(stockTimePlot, netflowTimePlot, availTimePlot, growthrateTimePlot, cols = 2)

```




### 4.3 Modeling Constraints—A Non-renewable Stock

Many systems modeled with system dynamics methods involve non-renewable stocks, which are stocks that when depleted cause another stock’s value to decline. The model below depicts how oil well growth, capital, is dependent on the availability of oil, the resource. As capital grows so does the amount of oil extracted. This process has several loops.



[Limits to growth for capital, constrained by a non-renewable resource](fig3.png)




There are five loops in the oil well model of varying types. The positive feedback loop R1 is driven by a desired investment in the capital, the oil wells, whose growth is determined by the growth fraction. Growth is measured by the number of oil wells; more wells equates to more extraction, capital, investment, and profit. The negative feedback loop, B1, accounts for the capital depreciation rate. Over time oil well equipment will need repairs or replacement thus driving depreciation. In the negative feedback loop B2, oil well equipment purchase (startup) costs are driven by the total revenue, amount of reinvestments, profit, and capital costs. More capital results in more oil extraction. However, more investments in capital reduces profits and lower available investments. The negative feedback loops B3 and B4 describe the final aspect of the model, the non-renewable resource oil. As extraction increases with capital, the resource eventually decreases along with extraction efficiency and rate. Efficiency is reduced due to the lack of natural pressure pushing the oil to the surface.

### 4.3.1 Equations

The capital stock has an initial value of 5 and integrates the difference of investments and depreciation (constant 5%).

$$Capital=\int _5^n Investments∗Depreciation$$

$$Depreciation=Capital∗DepreciationRate$$

    DepreciationRate=0.05


The amount of desired investment is the product of growth fraction and capital.

$$DesiredInvestment=DesiredGrowthFraction∗Capital$$

$$DesiredGrowthFraction=0.07$$


Growth is limited by the resource since efficiency decreases, which is defined by extraction rate and is important to understanding feedback loops B3 and B4. The plot of resource by efficiency (below) displays this relationship.

$$Resource=\int_n1000−Extraction$$

$$Extraction=Capital∗ExtractionEfficiencyPerUnitCapital$$


With capital costs set at 10% of capital and revenue per unit extracted at 3, the total revenue can be determined and then deducted to attain the profit.

$$TotalRevenue=RevenuePerUnitExtracted∗Extraction$$

$$RevenuePerUnitExtracted=3$$

$$CapitalCosts=Capital∗0.10$$

$$Profit=TotalRevenue−CapitalCosts$$


The maximum capital investment is calculated by allocating $12\%$ of profits reinvested at a cost of 2 per investment.

$$FractionProfitsReinvested=0.12$$

$$CapitalFunds=Profit∗FractionProfitsReinvested$$

$$CostPerInvestment=2$$

$$MaximumInvestment=CapitalFundsCostPerInvestment$$


The MIN function is applied as a practical verification for maximum possible investment value so as to not exceed a particular goal.

$$Investment=MIN(DesiredInvestment,MaximumInvestment)$$


The model is initialized with a time period and step along with stocks and auxiliary values and rates.

```{r}
# Set the time period and step
START <- 0
FINISH <- 200
STEP <- 0.25

# Set the stocks and auxiliaries
simtime <- seq(START, FINISH, by = STEP)
stocks <- c(sCapital = 5, sResource = 1000)
auxs <- c(aDesired.Growth = 0.07,
          aDepreciation = 0.05,
          aCost.Per.Investment = 2.00,
          aFraction.Reinvested = 0.12,
          aRevenue.Per.Unit = 3.00)
```



The **approxfun()** function interpolates input and output vectors (resource and extraction) in order to “linearize” the new function’s output, extraction efficiency per unit capital:

ExtractionEfficiencyPerUnitCapital=

Resource |	Efficiency
---------|------------
1000|1.0
900	|0.99
800	|0.98
700	|0.96
600	|0.92
500	|0.85
400	|0.75
300	|0.63
200	|0.45
100	|0.25
0   |	0


```{r}
# Set the x and y vectors for linearization by the `approxfun` command
x.Resource <- seq(0, 1000, by = 100)
y.Efficiency <- c(0, 0.25, 0.45, 0.63, 0.75, 0.85, 0.92, 0.96, 0.98, 0.99, 1.0)

func.Efficiency <- approxfun(x = x.Resource,
                             y = y.Efficiency,
                             method = "linear",
                             yleft = 0, yright = 1.0)

plot(x.Resource, y.Efficiency)
```




### 4.3.2 Function
Given the previous equations and functions, the model’s function is constructed and the ode function called to create a data frame of the output.

```{r}
# Create the function for the model and the list to return for output
model <- function(time, stocks, auxs){
  with(as.list(c(stocks, auxs)), {
    aExtr.Efficiency <- func.Efficiency(sResource)
    fExtraction <- aExtr.Efficiency * sCapital
    aTotal.Revenue <- aRevenue.Per.Unit * fExtraction
    aCapital.Costs <- sCapital * 0.10
    aProfit <- aTotal.Revenue - aCapital.Costs
    aCapital.Funds <- aFraction.Reinvested * aProfit
    aMaximum.Investment <- aCapital.Funds / aCost.Per.Investment
    aDesired.Investment <- sCapital * aDesired.Growth
    fInvestment <- min(aMaximum.Investment,
                       aDesired.Investment)
    fDepreciation <- sCapital * aDepreciation
    dS_dt <- fInvestment - fDepreciation
    dR_dt <- -fExtraction
    
    return(list(c(dS_dt, dR_dt),
                CapitalCosts = aCapital.Costs,
                Capital = sCapital,
                DesiredInvestment = aDesired.Investment,
                MaximumInvestment = aMaximum.Investment,
                Investment = fInvestment,
                Depreciation = fDepreciation,
                Extraction = fExtraction))
  })
}

```


```{r}
# Create the data frame to use for visualizations
nonNewModel <- data.frame(ode(y = stocks, times = simtime, func = model,
                    parms = auxs, method = "euler"))
```



### 4.3.3 Outputs

Inspecting the plot generated below depicting extraction, capital, resources, and depreciation shows the relationships of the stocks and flows. Exponential extraction rates precedes capital while reducing the resources, increasing investments (black dashed line) and capital depreciation. Ultimately the resource becomes too expensive to extract as does the capital overhead causing capital and profits to decline.

```{r}
# Plot the output using `multiplot`
extraction <- nonNewModel %>% 
  ggplot() +
  geom_line(aes(time, Extraction), color = "orange") +
  ylab("Extraction Rate") +
  xlab("Year")

capital <- nonNewModel %>% 
  ggplot() +
  geom_line(aes(time, sCapital), color = "red") +
  ylab("Capital") +
  xlab("Year")

resources <- nonNewModel %>% 
  ggplot() +
  geom_line(aes(time, sResource), color = "purple") +
  ylab("Resource") +
  xlab("Year")

depreciation <- nonNewModel %>% 
  ggplot() +
  geom_line(aes(time, Depreciation), color = "green") +
  geom_line(aes(time, Investment), linetype = "dashed", size = 0.75) +
  ylab("Depreciation") +
  xlab("Year")

multiplot(extraction, capital, resources, depreciation, cols = 2)

```





### 4.3.4 Maximum Values

For this scenario, the years during which maximum capital and extraction rates occur are calculated. These are particularly useful when conveying strengths and liabilities.

```{r}
# Find the time (year) in which the maximum capital and extraction rate are realized
nonNewModel[which.max(nonNewModel$sCapital), "time"]
## [1] 87
nonNewModel[which.max(nonNewModel$Extraction), "time"]
```





### 4.3.5 Exercise


Increase the capital depreciation to 25% and plot the output using multiplot. Explain the results.

```{r}
# Change the aDepreciation from 0.05 to 0.25
auxs <- c(aDesired.Growth = 0.07,
          aDepreciation = 0.25,
          aCost.Per.Investment = 2.00,
          aFraction.Reinvested = 0.12,
          aRevenue.Per.Unit = 3.00)

# Run the `model2` function
model2 <- function(time, stocks, auxs){
  with(as.list(c(stocks, auxs)), {
    aExtr.Efficiency <- func.Efficiency(sResource)
    fExtraction <- aExtr.Efficiency * sCapital
    aTotal.Revenue <- aRevenue.Per.Unit * fExtraction
    aCapital.Costs <- sCapital * 0.10
    aProfit <- aTotal.Revenue - aCapital.Costs
    aCapital.Funds <- aFraction.Reinvested * aProfit
    aMaximum.Investment <- aCapital.Funds / aCost.Per.Investment
    aDesired.Investment <- sCapital * aDesired.Growth
    fInvestment <- min(aMaximum.Investment,
                       aDesired.Investment)
    fDepreciation <- sCapital * aDepreciation
    dS_dt <- fInvestment - fDepreciation
    dR_dt <- -fExtraction
    
    return(list(c(dS_dt, dR_dt),
                CapitalCosts = aCapital.Costs,
                Capital = sCapital,
                DesiredInvestment = aDesired.Investment,
                MaximumInvestment = aMaximum.Investment,
                Investment = fInvestment,
                Depreciation = fDepreciation,
                Extraction = fExtraction))
  })
}

# Run the `data.frame` command
nonNewModel2 <- data.frame(ode(y = stocks, times = simtime, func = model2,
                    parms = auxs, method = "euler"))

# Run `ggplot` for each variable and `multiplot` for the visualization
extraction <- nonNewModel2 %>% 
  ggplot() +
  geom_line(aes(time, Extraction), color = "orange") +
  ylab("Extraction Rate") +
  xlab("Year")

capital <- nonNewModel2 %>% 
  ggplot() +
  geom_line(aes(time, sCapital), color = "red") +
  ylab("Capital") +
  xlab("Year")

resources <- nonNewModel2 %>% 
  ggplot() +
  geom_line(aes(time, sResource), color = "purple") +
  ylab("Resource") +
  xlab("Year")

depreciation <- nonNewModel2 %>% 
  ggplot() +
  geom_line(aes(time, Depreciation), color = "green") +
  geom_line(aes(time, Investment), linetype = "dashed", size = 0.75) +
  ylab("Depreciation") +
  xlab("Year")

multiplot(extraction, capital, resources, depreciation, cols = 2)

```



### 4.3.6 Creating Scenarios: Effect of Growth Rate on Results
System dynamics modeling often allows users to investigate multiple secarios for consideration. Using the rbind() function to append data sets, a range of growth rates is then displayed in the graph below. Given the outputs of the model (above), the changes to capital, extraction rate, and investments becomes apparent. Typically, the more aggressive the growth rate the more rapid the extraction and subsequently the consumption of the resource.


```{r}
# Change the aDepreciation back to 0.05
auxs <- c(aDesired.Growth = 0.07,
          aDepreciation = 0.05,
          aCost.Per.Investment = 2.00,
          aFraction.Reinvested = 0.12,
          aRevenue.Per.Unit = 3.00)

# Function to automate `rbind()` below. The function formals include `stocks`, `simtime`, `model`, and `auxs`, which were defined in the arguments of the previous data frame command. 

scenarios <- function(GRate, n, stocks, simtime, model, auxs, step = 0.01) {
    vec <- seq(GRate, by = step, length.out = n)
    base <- data.frame()
    for (i in 1:length(vec)) {
        auxs["aDesired.Growth"] <- vec[i]
        temp <- data.frame(ode(y = stocks, times = simtime, func = model, parms = auxs, method = "euler"))
        temp$GR <- paste0("GR = ",100*vec[i],"%")
        base <- rbind(base,temp)
    }
    return(base)
}

base <- scenarios(0.05, 7, stocks, simtime, model, auxs)

# Plot the scenarios
base %>% 
  ggplot() +
  geom_line(aes(time, Extraction, color = base$GR)) +
  xlab("Year") +
  ylab("Extraction Rate") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(title = NULL)) +
  ggtitle("Growth Rate Scenarios")


```



### 4.3.7 Exercise

Plot a set of six growth rate scenarios from two to seven percent. Identify the year the maximum extraction rate occurs.

```{r}
# Generate the growth rate scenarios
base <- scenarios(0.02, 6, stocks, simtime, model, auxs)

# Plot the scenarios
base %>% 
  ggplot() +
  geom_line(aes(time, Extraction, color = base$GR)) +
  xlab("Year") +
  ylab("Extraction Rate") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(title = NULL)) +
  ggtitle("Growth Rate Scenarios")





```

```{r}
# Determine the year maximum extraction rate occurs
base[which.max(base$Extraction), "time"]
```



It is important to acknowledge that this model is paramount in a variety of industries where growth has to be considered based on the available resources, whether that be oil, patients, or customers where several independent variables act upon a dependent variable. Given the toy instance described above, it is easy to imagine the complexity of a more realistic and pragmatic model with dozens of stocks and flows.

## References

Duggan, Jim. 2016. System Dynamics Modeling with R: Lecture Notes in Social Networks.

Forrester, Jay W. 1968. Principles of Systems. Cambridge: Wright-Allen Press, Inc.

System Dynamics Society. 2017. “Origin of System Dynamics.” Accessed August 2. http://www.systemdynamics.org/DL-IntroSysDyn/origin.htm.

