Fibonacci Rabbits

This is a small, discrete time model that explores the physical interpretation of the Fibonacci sequence. See my blog post about this model for details.

Fibonacci2.vpm This runs with Vensim PLE, but users might want to use the Model Reader in order to load the included .cin file with non-growing eigenvector settings.

Oscillation from a purely positive loop

Replicated by Mohammad Mojtahedzadeh from Alan Graham’s thesis, or created anew with the same inspiration. He created these models in the course of his thesis work on structural analysis through pathway participation matrices.

Alan Graham, 1977. Principles on the Relationship Between Structure and Behavior of Dynamic Systems. MIT Thesis. Page 76+

These models are pure positive feedback loops that don’t exhibit exponential growth (under the right initial conditions). See my blog post for a discussion of the details.

These are generic models, and therefore don’t have units. All should run with Vensim PLE, except the generic gain matrix version which uses arrays and therefore requires an advanced version or the Model Reader.

The original 4th order model, replicated from Alan’s thesis: PurePosOscill4.vpm – note that this includes a .cin file with an alternate stable initialization.

My slightly modified version, permitting initialization with different gains at each level: PurePosOscill4alt.vpm

Loops of different orders: 3.vpm 6.vpm 8.vpm 12.vpm (I haven’t spent much time with these. It appears that the high-order versions transition to growth rather quickly – my guess is that this is an artifact of numerical precision, i.e. any tiny imprecision in the initialization introduces a bit of the growth eigenvector, which quickly swamps the oscillatory signal. It would be interesting to try these in double precision Vensim to see if I’m right.)

Stable initializations: 2stab.vpm 12stab.vpm

A generic version, representing a system as a generic gain matrix, so you can use it to explore any linear unforced variant: Generic.vpm

Bifurcating Salmon

A nifty paper on nonlinear dynamics of salmon populations caught my eye on ArXiv.org today. The math is straightforward and elegant, so I replicated the model in Vensim.

A three-species model explaining cyclic dominance of pacific salmon

Authors: Christian Guill, Barbara Drossel, Wolfram Just, Eddy Carmack

Abstract: The four-year oscillations of the number of spawning sockeye salmon (Oncorhynchus nerka) that return to their native stream within the Fraser River basin in Canada are a striking example of population oscillations. The period of the oscillation corresponds to the dominant generation time of these fish. Various – not fully convincing – explanations for these oscillations have been proposed, including stochastic influences, depensatory fishing, or genetic effects. Here, we show that the oscillations can be explained as a stable dynamical attractor of the population dynamics, resulting from a strong resonance near a Neimark Sacker bifurcation. This explains not only the long-term persistence of these oscillations, but also reproduces correctly the empirical sequence of salmon abundance within one period of the oscillations. Furthermore, it explains the observation that these oscillations occur only in sockeye stocks originating from large oligotrophic lakes, and that they are usually not observed in salmon species that have a longer generation time.

The paper does a nice job of connecting behavior to structure, and of relating the emergence of oscillations to eigenvalues in the linearized system.

Units balance, though I had to add a couple implicit scale factors to do so.

The general results are qualitatitively replicable. I haven’t tried to precisely reproduce the authors’ bifurcation diagram and other experiments, in part because I couldn’t find a precise specification of numerical methods used (time step, integration method), so I wouldn’t expect to succeed.

Unlike most SD models, this is a hybrid discrete-continuous system. Salmon, predator and zooplankton populations evolve continuously during a growing season, but with discrete transitions between seasons.

The model uses SAMPLE IF TRUE, so you need an advanced version of Vensim to run it, or the free Model Reader. (It should be possible to replace the SAMPLE IF TRUE if an enterprising person wanted a PLE version). It would also be a good candidate for an application of SHIFT IF TRUE if someone wanted to experiment with the cohort age structure.

sockeye.vmf

For a more policy-oriented take on salmon, check out Andy Ford’s work on smolt migration.

DICE

This is a replication of William Nordhaus’ original DICE model, as described in Managing the Global Commons and a 1992 Science article and Cowles Foundation working paper that preceded it.

There are many good things about this model, but also some bad. If you are thinking of using it as a platform for expansion, read my dissertation first.

Units balance.

I provide several versions:

  1. Model with simple heuristics replacing the time-vector decisions in the original; runs in Vensim PLE
  2. Full model, with decisions implemented as vectors of points over time; requires Vensim Pro or DSS
  3. Same as #2, but with VECTOR LOOKUP replaced with VECTOR ELM MAP; supports earlier versions of Pro or DSS
    • DICE-vec-6-elm.mdl (you’ll also want a copy of DICE-vec-6.vpm above, so that you can extract the supporting optimization control files)

Note that there may be minor variances from the published versions, e.g. that transversality coefficients for the state variables (i.e. terminal values of the states for optimization) are not included. The optimizations use fewer time decision points than the original GAMS equivalents. These do not have any significant effect on the outcome.

Java Vensim helper

MIT’s Climate Collaboratorium has posted java code that it used to wrap C-LEARN as a web service using the multicontext .dll. If you’re doing something similar, you may find the code useful, particularly the VensimHelper class. The liberal MIT license applies. However, be aware that you’ll need a license for the Vensim multicontext .dll to go with it.

Cumulative Normal Distribution

Vensim doesn’t have a function for the cumulative normal distribution, but it’s easy to implement via a macro. I used to use a polynomial cited in Numerical Recipes (error function, Ch. 6.2):

:MACRO: NCDF(x)
NCDF = 1-Complementary Normal CDF
~	dmnl
~		|
Complementary Normal CDF=
ERFCy/2
~	dmnl
~		|
ERFCy = IF THEN ELSE(y>=0,ans,2-ans)
~	dmnl
~	http://www.library.cornell.edu/nr/bookcpdf/c6-2.pdf
|
y = x/sqrt(2)
~	dmnl
~		|
ans=t*exp(-z*z-1.26551+t*(1.00002+t*(0.374092+t*(0.0967842+
t*(-0.186288+t*(0.278868+t*(-1.1352+t*(1.48852+
t*(-0.822152+t*0.170873)))))))))
~	dmnl
~		|
t=1/(1+0.5*z)
~	dmnl
~		|
z = ABS(y)
~	dmnl
~		|
:END OF MACRO:

I recently discovered a better approximation here, from algorithm 26.2.17 in Abromowitz and Stegun, Handbook of Mathematical Functions:

:MACRO: NCDF2(x)
NCDF2 =  IF THEN ELSE(x >= 0,
(1 - c * exp( -x * x / 2 ) * t *
( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 )),  ( c * exp( -x * x / 2 ) * t *
( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 ))
)
~     dmnl
~     From http://www.sitmo.com/doc/Calculating_the_Cumulative_Normal_Distribution
Implements algorithm 26.2.17 from Abromowitz and Stegun, Handbook of Mathematical 
Functions. It has a maximum absolute error of 7.5e^-8.
http://www.math.sfu.ca/
|
c  =  0.398942
~     dmnl
~           |
t = IF THEN ELSE( x >= 0, 1/(1+p*x), 1/(1-p*x))
~     dmnl
~           |
b5 =  1.33027
~     dmnl
~           |
b4 = -1.82126
~     dmnl
~           |
b3 =  1.78148
~     dmnl
~           |
b2 = -0.356564
~     dmnl
~           |
b1 =  0.319382
~     dmnl
~           |
p  =  0.231642
~     dmnl
~           |
:END OF MACRO:

In advanced Vensim versions, paste the macro into the header of your model (View>As Text). Otherwise, you can implement the equations inside the macro directly in your model.

Rental car stochastic dynamics

This is a little experimental model that I developed to investigate stochastic allocation of rental cars, in response to a Vensim forum question.

There’s a single fleet of rental cars distributed around 50 cities, connected by a random distance matrix (probably not physically realizable on a 2D manifold, but good enough for test purposes). In each city, customers arrive at random, rent a car if available, and return it locally or in another city. Along the way, the dawdle a bit, so returns are essentially a 2nd order delay of rentals: a combination of transit time and idle time.

The two interesting features here are:

  • Proper use of Poisson arrivals within each time step, so that car flows are dimensionally consistent and preserve the integer constraint (no fractional cars)
  • Use of Vensim’s ALLOC_P/MARKETP functions to constrain rentals when car availability is low. The usual approach, setting actual = MIN(desired, available/TIME STEP), doesn’t work because available is subscripted by 50 cities, while desired has 50 x 50 origin-destination pairs. Therefore the constrained allocation could result in fractional cars. The alternative approach is to set up a randomized first-come, first-served queue, so that any shortfall preserves the integer constraint.

The interesting experiment with this model is to lower the fleet until it becomes a constraint (at around 10,000 cars).

Documentation is sparse, but units balance.

Requires an advanced Vensim version (for arrays) or the free Model Reader.

carRental.vpm carRental.vmf

Update, with improved distribution choice and smaller array dimensions for convenience:

carRental2.mdl carRental2.vpm

Pink Noise

In a continuous time dynamic model, representing noise as a random draw at every time step can be problematic. As the time step is decreased, the high frequency power of the noise spectrum increases accordingly, potentially changing the behavior. In the limit of small time steps, the resulting white noise has infinite power, which is not physically realistic.

The solution is to use pink noise, which is essentially white noise filtered to cut off high frequencies. SD models from the bad old days typically employed a pink noise generating structure that employed uniformly distributed white noise, relying on the central limit theorem to yield a normally distributed output. Ed Anderson improved that structure to incorporate a normally distributed input, which works better, especially if the cutoff frequency is close to the inverse of the time step.

Two versions of the model are attached: one for advanced versions of Vensim, which permit implementation as a :MACRO:, for efficient reuse. The other works with Vensim PLE.

PinkNoise2010.mdl PinkNoise2010.vmf PinkNoise2010.vpm

PinkNoise2010-PLE.vmf PinkNoise2010-PLE.vpm

Contributed by Ed Anderson, updated by Tom Fiddaman

Notes (also in the model files):

Description: The pink noise molecule described generates a simple random series with autocorrelation. This is useful in representing time series, like rainfall from day to day, in which today’s value has some correlation with what happened yesterday. This particular formulation will also have properties such as standard deviation and mean that are insensitive both to the time step and the correlation (smoothing) time. Finally, the output as a whole and the difference in values between any two days is guaranteed to be Gaussian (normal) in distribution.

Behavior: Pink noise series will have both a historical and a random component during each period. The relative “trend-to-noise” ratio is controlled by the length of the correlation time. As the correlation time approaches zero, the pink noise output will become more independent of its historical value and more “noisy.” On the other hand, as the correlation time approaches infinity, the pink noise output will approximate a continuous time random walk or Brownian motion. Displayed above are two time series with correlation times of 1 and 8 months. While both series have approximately the same standard deviation, the 1-month correlation time series is less smooth from period to period than the 8-month series, which is characterized by “sustained” swings in a given direction. Note that this behavior will be independent of the time-step. The “pink” in pink noise refers to the power spectrum of the output. A time series in which each period’s observation is independent of the past is characterized by a flat or “white” power spectrum. Smoothing a time series attenuates the higher or “bluer” frequencies of the power spectrum, leaving the lower or “redder” frequencies relatively stronger in the output.

Caveats: This assumes the use of Euler integration with a time step of no more than 1/4 of the correlation time. Very long correlation times should be avoided also as the multiplication in the scaled white noise will become progressively less accurate.

Technical Notes: This particular form of pink noise is superior to that of Britting presented in Richardson and Pugh (1981) because the Gaussian (Normal) distribution of the output does not depend on the Central Limit Theorem. (Dynamo did not have a Gaussian random number generator and hence R&P had to invoke the CLM to get a normal distribution.) Rather, this molecule’s normal output is a result of the observations being a sum of Gaussian draws. Hence, the series over short intervals should better approximate normality than the macro in R&P.

MEAN: This is the desired mean for the pink noise.

STD DEVIATION: This is the desired standard deviation for the pink noise.

CORRELATION TIME: This is the smooth time for the noise, or for the more technically minded this is the inverse of the filter’s cut-off frequency in radians.