EGR 312 Final Project Example

Problem Description

The problem described below was adapted from a SIMIODE modeling challenge. SIMIODE is a community that provides instructor resources for teaching differential equations through modeling. Parts of the original problem description were taken from an existing SIMIODE resource.

In 2011 the country of Sudan split into Sudan and South Sudan. The split followed many years of strife. Unfortunately, the split itself led to further violence, which in return resulted in the creation of refugee settlements in Uganda. There are many aspects of the dynamics of the lives of the people within these camps that are not understood, and the camp settlements themselves are the subject of observation and research.

The settlements in Uganda were originally formed by different outside organizations. Each organization tries to impose its own methods to resolve disputes that may arise between people in the settlement. For example, some organizations may set up a system in which people are expected to negotiate with one another, while another organization may prefer to create a system in which outside arbiters make decisions by which others must abide. Another common approach is to require the use of a mediator who acts as a guide to help the parties involved find their own solution to which all can agree.

In a study of some of the refugees who fled South Sudan and moved to settlements in Uganda, it was found that many people also brought their own systems and preferences with them and sought to incorporate their own traditions into the process. I can use a mathematical model of this system to simulate different scenarios, with the goal of informing policy and management strategies. Questions to explore include the following: What happens when multiple practices and traditions are brought together? How do the different practices come together and create a new system in the context of a refugee settlement? How does a process of dispute resolution change and evolve?

Model Conceptualization

I can simplify this complex system by focusing on the three methods described above: negotiation, arbitration, and mediation and make some assumptions:

1. All three dispute resolution strategies are available in a settlement.

2. Each person initially has a specific preference for mediation, negotiation, and arbitration.

3. When disputes occur, they are settled using only one of the three methods.

4. When a dispute occurs between two people with different initial preferred strategies, the method that is chosen to resolve the dispute becomes the preferred strategy for both people.

5. Disputes occur every time people with different initial preferences interact.

6. Only a fraction of the population chooses to interact with those that have different preferences.

7. Each settlement has a static population.

A visual depiction of the simplified system under these assumptions is provided in Figure 2. This schematic supported the development of the mathematical model presented below.

Explanation and Justification of Approaches

In Part 1 of this project, I simulated the model using both a simple and complex numerical method for simulating differential equations. The simple approach was Euler's method, a 1^st order forward finite difference approximation of the differential terms. While this is the simplest method, it also has the largest truncation error. I compared results from Euler's method with the 5th-order Runge Kutta method, which has lower truncation error but introduces more round-off error due to the additional calculations that need to be implemented. Prior to implementing these simulations, I had to identify values for each of the four model parameters. I came up with three alternative scenarios that I thought might occur in a refugee camp, with the goal of producing realistic simulations.

Scenario #1:

I will assume that 20% of the population chooses to interact with people who have different conflict resolution strategies than their own. This assumption is based on my own qualitative perception of how my peers tend to spend about 80% of their time within their own social circles but branch out and interact with others about 20% of the time. For the fractions of disputes that are resolved by each group, I assumed that arbitrators would be least comfortable with negotiators and would typically insist on bringing in an arbitrator. However, I assumed that arbitrators would be very open to mediation. I believe mediators would more often insist on a mediator being present but would also acquiesce to the negotiation strategy in many scenarios. Based on these assumptions, the parameter values for scenario #1 are defined as follows:

k = 0.2

aNA = 0.7 (nNA = 0.3)

nMN = 0.4 (mMN = 0.6)

mMA = 0.8 (aMA = 0.2)

Scenario #2

For scenario #2, I assumed that the three groups are more inclined to stay within their own groups and avoid interacting with others unless necessary. Because these groups are avoiding interactions with one another, I assumed that they would want to resolve disputes quickly and would be more receptive to a negotiation strategy to potentially handle conflict in an efficient manner. Similarly, arbitration would be preferred over mediation to ensure that a decision is made following the dispute. Based on these assumptions, the parameter values for scenario #3 are as follows:

k = 0.05

aNA = 0.4 (nNA = 0.6)

nMN = 0.8 (mMN = 0.2)

mMA = 0.3 (aMA = 0.7)

Scenario #3

For scenario #3, I assumed the groups interact more and have had more experience with conflict. I assumed the groups had many negative experiences following negotiation because they felt they were taken advantage of. They might also express discontent with the mediation strategy if the two groups do not resolve their conflict following the dispute. Therefore, I assumed the groups tend to prefer arbitration to bring in an outside party and ensure that a final decision is made.

k = 0.3

aNA = 0.75 (nNA = 0.25)

nMN = 0.4 (mMN = 0.6)

mMA = 0.4 (aMA = 0.6)

In Part 2, I developed an alternative model that did not have a static population (assumption #7). I used (simulated) data that represented the incoming population of each group to calibrate regression curves, using least squares regression. The resulting fitted equations were used as inputs into the system of differential equations to simulate the evolution of arbitrators, negotiators, and mediators in this revised system, shown in Part 3.

In Part 4, I used (simulated) data that represented the costs of running this refugee camp to estimate operational costs over time. To compare the total operating costs associated with different scenarios, I used numerical integration, where the estimated monthly costs are weighted by the populations of negotiators, arbitrators, and mediators and summed over time. This relatively simply approach is adequate considering the uncertainty ingrained in the model structure, population growth trends, and operating costs; additional errors introduced by the numerical methods are relatively minor.

Part 1: Model Simulation

Numerical Approach

Euler's Method

5th Order Runge Kutta

Pseudocode

• Define constant parameter values

• Define parameter values that vary for each scenario

• Set initial conditions

• Simulate system using Euler's method

  • for each time step

    • simulate N, A, and M population

  • end loop

• Calculate truncation error from Euler's method for each time step

• Simulate system using 5th order RK

• for each time step

  • calculate k1 for N, A, and M

  • calculate k2 = f(k1) for N, A, and M

  • calculate k3 = f(k1, k2) for N, A, and M

  • calculate k4 = f(k2, k3) for N, A, and M

  • calculate k5 = f(k1, k4) for N, A, and M

  • calculate k6 = f(k1, k2, k3, k4, k5) for N, A, and M

  • simulate system using k1,k2,k3,k4,k5,k6 for N, A, and M

• end loop

• plot results

  • Euler's and 5th order RK together

  • truncation error

  • difference between methods compared to truncation error

Results

This simulation was applied for each of the three scenarios identified and the changes in population of each group over time are illustrated in Figure 3 for Euler’s method and Figure 4 for the 5^th order Runge-Kutta (Butcher’s) method. The two numerical approaches are compared for Scenario #1 in Figure 5.

Error Analysis

The types of error associated with our numerical methods include truncation error and round-off error. The truncation error is associated with the approximation of the derivative and should be larger for Euler’s method compared to Butcher’s method. We can use Taylor series expansion to calculate the truncation error associated with Euler’s method, which is a first order forward finite difference approximation and is equivalent a first order Taylor Series expansion. Therefore, the higher order terms in the Taylor series expansion provide the truncation error from Euler’s method. In our system, the third derivative of N(t), M(t) and A(t) is zero, so only the second derivative terms needs to be computed to quantify truncation error. A sample calculation for computing the truncation error associated with estimating N(1) from N(0) in Scenario #1 is provided below. The truncation error for the simulation of N, M, and A for time steps for scenario #1 was computed and plotted in Figure 6. This error is not significant, especially considering the magnitude is in fractions of people, which would not exist in reality.

Round-off error can be quantified, although it would be difficult as it would require tracking the number of flops and considering how the computer stores numbers and its memory limitations. However, we know that Butcher’s method should have lower truncation error but required more computations so we can examine the difference between simulations from Euler’s and Butcher’s methods to view round-off error qualitatively. The differences between the two simulations are plotted for Scenario #1 in Figure 7 and compared to the computed truncation error. Note that the difference between the two methods is significantly greater than the truncation error. This indicates that round-off error dominates and the added computations from Butcher’s are not justified. Euler’s method should perform well based on the plot of the shape of the functions in Figure 3, since they are nearly linear for individual time steps. Another strategy for reducing truncation error would be to reduce the step size, but this would also increase round-off error.

Discussion

From this initial model simulation, we can see that one dispute resolution strategy will eventually become the dominant strategy in the refugee camp. This observation is logical considering the mathematical model; as long as the differences in the fraction parameters are greater than zero, the population with the most preferred strategy will continue to grow. The k parameters will slow down or speed up the rate at which groups adopt alternative resolution strategies. Therefore, our assumptions on the parameter values have a very substantial impact on the simulation and determine which group will become the dominant population. The differences between the numerical simulations were fairly significant, especially considering the truncation error was insignificant. However, considering the assumptions made in this mathematical model, this simulation is not very realistic and there are large uncertainties that we cannot quantify without real data on how these populations actually evolve over time. Considering the larger uncertainties of the model, and the observation that truncation error is not very significant, Euler’s method is the most practical approach to simulate this simplistic model.

Part 2: Alternative Model with Population Growth

If we wanted to re-examine the last assumption (#7: Each settlement has a static population) and had some data on the increase in the population of each settlement, we could adjust the system of equations according to Equations 4 - 6, illustrated in Figure 8.

Numerical Approach

I was provided with data on the incoming populations to each community. Using curve fitting techniques, functions can be fit to this data and the more complex model described in Equations 4, 5, and 6.

Before fitting a curve, each dataset should be plotted to identify appropriate regression models. Scatter plots for each group are shown in Figure 7. The group pattern for the Negotiators appears linear, the mediator growth appears to follow a saturation growth curve, and the arbitrator growth resembles a sinusoid. Each of these functions were fit to the data using the generalized linear least squares approach, illustrated next.

In the saturation growth rate model, note that the coefficients need to be converted back to the nonlinear form after solving for a0 and a1.

To confirm that this data follows the saturation growth model, I plotted the 1/y versus 1/x to check that there was a linear relationship once the model is inverted. From this plot, there was a clear outlier. This point was therefore ignored when fitting the model.

Pseudocode

• plot scatter plots of each variable to identify appropriate models

• fit linear model to negotiators

  • set up Z matrix with intercept and slope term

  • solve for coefficients using inv(Z'Z)(Z'Y)

  • get fitted model using coefficients

• calculate total sum of squares, sum of squared residuals, R-squared, and standard error of Y|X

• plot data versus model

• fit saturation growth rate model

  • set up Z matrix using intercept and 1/X

  • solve for coefficients using inv(Z'Z)(Z'Y)

  • convert fitted coefficients to model coefficients

  • get fitted model using coefficients

• calculate total sum of squares, sum of squared residuals, R-squared, and standard error of Y|X

• plot data versus model

• fit sinusoid model

  • set up Z matrix using intercept, cos(omega*x), sin(omega*x)

  • solve for coefficients using inv(Z'Z)(Z'Y)

  • convert fitted coefficients to model coefficients - a, c, theta

  • get fitted model using coefficients

• calculate total sum of squares, sum of squared residuals, R-squared, and standard error of Y|X

• plot data versus model

Results

Error Analysis

Key error metrics associated with these models are the R-squared values and the standard error. The R-squared represents how much of the variation, or spread, in the population growth data is explained by the fitted model. An R-squared value close to 1 would indicate that the fit was very good and the population growth trends over time can be accurately predicted (assuming stationarity). The standard error is calculated from the sum of squared residuals and provides an estimate of the error in population growth for an instance in time. Assuming our errors are normally distributed, one standard error away from our best guess (the fitted model) provides the 68% confidence interval. This interval is plotted in Figures 8, 9, and 10 in the results section for each model fit. The standard error and r-squared values were computed according to the equations below, where y_bar represents the mean of the observations and y_hat is the best guess of y_i using the fitted model. The computed values are summarized in Table 1.

Discussion

Based on the errors reported in Table 1, all of the models are acceptable and appear to fit the datasets well. The Mediation curve has the lowest error and highest R-squared, and the simple linear model has the most error and lowest R-squared. The fitted models can be used to extrapolate the population growth past the period of available data (100 months) assuming the relationship stays the same in the future. This is helpful because it allows us to perform the model simulation over longer time periods. The trends of all three populations is shown in Figure 13 over a duration of 100 months. Note that the negotiator group will continue to grow if the data is extrapolate, which is not realistic based on the capacity of the refugee camp. These relationships would likely change over time. We could also consider people dying and leaving the camp to improve the model.

Part 3: Revised Simulation with Population Growth

A sample calculation for computing N(1), M(1), and A(1) from N(0), M(0), and A(0) for scenario #1 is shown below, where the fitted equations from part 2 have been applied to NIN(t), MIN(t), and AIN(t). Note that in the next time step, Ptot(t) will change and will no longer be constant at 1500.

Pseudocode

• create time series of incoming populations using fitted regression models

• plot results to conceptualize and logic check

• create total population time series

• repeat Euler's simulation using incoming population time series

• plot population over time for each scenario

• evaluate uncertainty from regression models

  • add +/- 1 std dev to each incoming population to capture 68% confidence interval

  • run high and low ends through simulation

  • plot results, showing upper and lower bounds

Results

The model was then simulated for all three scenarios and the results are shown in Figure 13. Note that scenario #2 is simulated over a longer time period than scenario #1 and #3 in order to capture more variation over time given the smaller k value.

Error Analysis

Truncation error and round-off error are present in this simulation, similar to the original model. However, these errors seem insignificant not when compared to the differences that result from changing the model structure (adding in population growth). In this simulation, I added an additional source of uncertainty – the population growth data. Each of the models we fit had some standard error associated with them. To evaluate this uncertainty, the model was simulated with the upper and lower bound of the 68% confidence interval of the population data. The resulting simulations for a period of 300 months are shown in Figure 15. Visually it is clear that the uncertainty from curve fitting is more substantial than the truncation error introduced through Euler’s method, which was plotted in part 1.

Discussion

In the previous simulation (refer to Figure 3) one of the three groups became the dominant population and approached the total population value, 1500. Now the populations can continue to increase. Despite the differences in incoming populations, the parameter values chosen for the fractions of disputes that are selected for a particular conflict resolution strategy appear to dictate which of the three strategies dominate long term. The other populations will not go to zero because the incoming population is positive each month, but these new populations will be converted to another strategy. If the population is growing, including this information in the model is very important to more accurately estimate how the populations will grow over time. However, I believe it is misleading to extrapolate the population growth data out to far in time, considering there are many external factors that would likely alter the temporal trends. Additionally, the refugee camp cannot continue to grow larger over time given physical and financial constraints associated with operating these camps.

Reflecting on the different assumptions made to develop these simulations, the structure of the conceptual model and the assumed parameters values have the greatest impact on the model and the conclusions I came to. The truncation and round-off error introduced by the numerical simulation of this system seems relatively insignificant for this scenario. However, if we were modeling a more exact system with less uncertainty, the error associated with numerical methods would be a more important consideration.

Part 4: Cost Analysis for Decision Making

The refugee camp managers must ensure they have enough funding to support each person in the camp. The cost to support each person varies monthly, as demands for water and electricity vary seasonally. The table below provides the monthly cost estimates. The managers would like to know the total operating cost for each scenario that I simulated.

The monthly operating costs can be computed over time by simply multiplying the monthly cost per person time the total population of the refugee camp. A vector of monthly costs was created in MATLAB by repeating the data in the table every 12 months.

To find the total cost over the 100 year period, numerical integration can be applied, where C(t) refers to the cost time series and Ptot(t) refers to the total population time series. Because we have functions for the population, a variety of numerical integration techniques could be applied to solve this integral, where higher accuracy is obtained for smaller interval sizes and successive applications of approximations. However, the monthly costs are provided as data so applying advanced integration methods would require interpolation between months or assumptions on the linearity of the function between data points. Therefore, the trapezoid method is an appropriate approach for finding the total cost. This approach is described below and accompanied by a sample calculation.

Costs associated with injuries and property damage tend to be higher for some groups compared to others, according to Table 3. If these costs are to be added to the total population costs, the populations associated with each scenario need to be considered. The monthly costs for each scenario were computed by simply multiply the monthly damage and injury costs by the populations for each group.

Numerical Approach

The refugee camp managers must ensure they have enough funding to support each person in the camp. The cost to support each person varies monthly, as demands for water and electricity vary seasonally. The table below provides the monthly cost estimates. The managers would like to know the total operating cost for each scenario that I simulated.

Pseudocode

• Input monthly cost data

• Calculate monthly operating costs for total population

• Plot results, showing costs from each population

• Calculate total cost using trapezoid rule

• Input damage and injury cost for each population

• Calculate cost from injury and damage for each population

• Calculate total cost over full time period using trapezoid rule

• Plot and compare costs.

Results

The monthly costs are plotted in Figure 16 for scenario #1. Note that the total cost is independent of the three scenarios since the total population is the same for all three. The seasonality of the costs can be observed in the plot, along with the upward trend caused by the population growth. This analysis demonstrates that the growing population in the camp is not sustainable over time.

The total cost over 100 months (over 4 years) was computed in MATLAB: Ctot = $143,500,000

Figure 17 shows the results for scenario #1. Despite the mediators being the dominant group, damage costs are highest from the negotiators. However, operating costs are still significantly higher.

The cost associated with damage, injuries and operating costs for each scenario are summarized in Table 4:

Discussion

According to the data on monthly expenses provided, scenario #2 would be the most expensive scenario, and scenario #3 would be the least expensive scenario; however, the differences between scenarios are relatively small considering the overall operating costs. Scenario #2 is the most expensive because the highest population of negotiators occurs, which tends to cause the most damage and injuries. Based on this information, the camp managers might put practices in place to reduce the population of negotiators. For instance, they could set up a partnership with an organization that finds arbitrators to serve refugee camps. If there are additional costs associated with these practices, this added cost would need to be included in the optimization problem.

In this case we are minimizing costs, which can be quantified according to the data provided along with the populations of each group. The constraints on the objective function would be the mathematical model simulation, which determines how the populations evolve and grow over time. Some additional trade-offs that might be added to the objective function would be measures of health, safety, and happiness. Perhaps a mental health survey could be distributed to these communities and poor indications of metal health could be translated to a cost that can be compared with the other financial costs. A challenge of this would be to quantify the value on mental health. Considering the very high operating costs, it is likely that this mental health cost would be very small compared to other financial expenses and would not play a major role in decision-making.

Despite the uncertainties and limitations of this modeling exercise, I did learn quite a bit about how these populations might evolve, and the numerical simulations triggered some important questions that refugee camps should consider. Numerical modeling gives us the opportunity to use simulate models that represent complex systems. In this exercise I looked at various scenarios that have different financial and social implications.