I introduced a recent paper, in the previous post, on the Tragedy of the Commons. The paper is published in Sustainability, an open access paper and is therefore available to everyone. One bit that I did not include in the paper, however, is the code that I used for the simulations therein. On the suggestion of a friend, I’ll publish them here in a series of posts. (Thanks Mauro)

The simulations are very simple, being based on difference equations, and were all coded in Octave. The code is not pretty, but I cleaned it up and added some comments. If you have questions, just drop me a comment on the blog. This first simulation is of a basic tragedy involving two users, and corresponds to Figures 1B-C in the paper. A resource is simulated to a steady state based on a Ricker model, for 100 steps, and then users begin to utilize it. They increase their utilization per time step at a steady rate (acceleration of utilization is constant), and the simulation is run until the resource is exhausted. The basic output of resource, one user’s benefit, and total utilization, are plotted in the figure here. The equations simulated are outlined in the excerpt from pg. 754 of the paper:
We now introduce a term to represent consumption or utilization by our human TOC agents.
\begin{aligned} R(t) &=& R(t)e^{r\left ( 1-\frac{R(t)}{K} \right )} - U(t)R(t) \nonumber \\ & \Rightarrow & R(t)\left [ e^{r\left ( 1-\frac{R(t)}{K} \right )} - U(t)\right ] \end{aligned}
where U(t) is the total fraction of R utilized by human users at time t, and ranges from 0 to 1. The term in square brackets on the right hand side of the equation is the modified growth rate of R. The resource is stable when this term is positive or zero, that is, the rate of resource renewal exceeds or is equal to utilization. Since U is the total standardized utilization rate of all users, it may be expanded to
$0 \leq U(t) = \sum_{i=1}^{N}u_{i}(t) \leq 1$
where there are N users and $u_{i}(t)$ is the standardized utilization rate of the $i^{th}$ user at time t. The benefit gained from utilization is modelled as
\begin{aligned} b_{i}(t) &=& f_{i}u_{i}(t)R(t)\nonumber \\ B(t) &=& R(t)\sum_{i=1}^{N}f_{i}u_{i} \end{aligned}
where $b_{i}$ is the benefit to user i, f is a factor that converts resources utilized into another commodity, for example converting harvested food to energy or currency, and B is the total benefit to all users.

The code follows, and can be executed by running it in an interactive Octave session. Please be aware that WordPress has somewhat limited options for formatting code. I’ve used Matlab formatting here, but wrapped lines might not be obvious. Single lines are all terminated with a semicolon, “;”, so just look out for those.

#RESOURCE DATA
#initial resource level
R0 = 1;
R = zeros(1,2);
R(1, 1) = 1;
R(1, 2) = 1;
#carrying capacity
K = 50;
#intrinsic growth rate r
r = 1.5;

#USER DATA
#There are 2 users
#initial acceleration of utilization
du0 = 1.05;
#individual user parameters, users 1 and 2
#initialize arrays
u1 = zeros(1,2);
u2 = zeros(1,2);
b1 = zeros(1,2);
b2 = zeros(1,2);
u1(1,1) = 1;
#initial utilization rate
u1(1,2) = .001;
#conversion factor of resource to benefit
f1 = 1;
b1(1,1) = 1;
b1(1,2) = 0;
u2(1,1) = 1;
#initial utilization rate
u2(1,2) = .001;
#conversion factor of resource to benefit
f2 = 1;
b2(1,1) = 1;
b2(1,2) = 0;
U = zeros(1,2);
U(1,1) = 1;
#total utilization
U(1,2) = b1(1,2) + b2(1,2);
du = zeros(1,2);
du(1,1) = 0;
du(1,2) = du0;

#BEGIN SIMULATION

#RUN RESOURCE TO STEADY STATE, NO UTILIZATION
for count2 = 2:100
R(count2, 1) = count2;
R(count2, 2) = R(count2-1, 2) * (exp(r*(1-(R(count2-1, 2)/K))));
u1(count2, 1) = count2;
u1(count2, 2) = u1(1,2);
u2(count2, 1) = count2;
u2(count2, 2) = u2(1,2);
du(count2,1) = count2;
du(count2,2) = du0;
endfor

#INITIATE UTILIZATION
for count1 = 101:265
R(count1, 1) = count1;
#modified Ricker model with utilization
R(count1, 2) = R(count1-1, 2) * (exp(r*(1-(R(count1-1, 2)/K))) - (u1(count1-1,2)+u2(count1-1,2)));
u1(count1, 1) = count1;
#increase utilization
u1(count1, 2) = u1(count1-1, 2) * du(count1-1,2);
#benefits (equal for both users in this case)
b1(count1, 1) = count1;
b1(count1, 2) = f1 * u1(count1-1, 2) * R(count1-1, 2);
b2(count1, 1) = count1;
b2(count1, 2) = f2 * u2(count1-1, 2) * R(count1-1, 2);
u2(count1, 1) = count1;
u2(count1, 2) = u2(count1-1, 2) * du(count1-1,2);
#total benefits
U(count1, 1) = count1;
U(count1, 2) = b1(count1, 2) + b2(count1, 2);
du(count1,1) = count1;
du(count1,2) = du(count1-1,2);
endfor