**EGR 312 Final Project Example**

**Problem Description**

**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?

Figure 1: Image of Sudenese refugee camp. Photo Credit: How refugees resolve disputes: insights from a Ugandan settlement. Vancluysen, Sarah and Bert Ingelaere. https://theconversation.com/how-refugees-resolve-disputes-insights-from-a-ugandansettlement-142508

**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.

Figure 2: Model schematic

**Explanation and Justification of Approaches**

**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**

**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.

Figure 3: Simulation of model #1 using Eulerâ€™s method for each scenario.

Figure 4: Simulation of Model #1 using Butcherâ€™s (4^{th} order Runge-Kutta method)

Figure 5: Comparison between Eulerâ€™s and Butcherâ€™s Methods for Scenario #1

## 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.

Figure 6: Truncation error calculated for Eulerâ€™s method

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.

Figure 7: Differences between alternative numerical simulations, compared to truncation 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**

**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.

Figure 8: Illustration of mathematical model with population growth

## 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.

Figure 9: Scatter plots of raw data showing population growth

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

Figure 10: Fitted linear model for negotiation population growth

Figure 11: Fitted saturation growth model for mediator population growth

Figure 12: Fitted sinusoid model for arbitrator population growth

## 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.

Figure 13: Population increases over time

**Part ****3****: ****Revised Simulation with Population Growth**

**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.

Figure 14: Simulation of Model #2 using Eulerâ€™s method for all three scenarios

## 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.

Figure 15: Propagation of error from curve fitting

## 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**

**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.

Figure 16: Monthly operating costs

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.

Figure 17: Scenario-based costs

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.

## MATLAB Code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% Prof Di Vittorio

% The purpose of this script is develop and test the final exam solution

% for EGR 312 - the dispute conflict resolution in South Sudan.

% Drafted March 16th 2021

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% model parameters

close all

clear all

%parameter definitions

% N = # of negotiators (DRC & Somalia)

% A = # of arbitrators (Burundi & Rwanda)

% M = # of mediators (South Sudan)

%% Scenario #1

sc = 1; %scenario

%proportion of disputes between N & A when N is chosen as resoltuion

%strategy

aNA = 0.1;

%proportion of disputes between N & A when A is chosen

nNA = 1- aNA;

%proportion of disputes between N & M when N is chosen

nNM = 0.2;

%propotion of disputes between N & M when M is chosen

mNM = 1 - nNM;

%propotion of disputes between M & A when M is chosen

mMA = 0.15;

%propotion of disputes between M & A when A is chosen

aMA = 1 - mMA;

%proportion of population that has disputes when randomly mixed

k = 0.05;

%time step

delt = 1; %months

T = 100; %number of months to simulate

%% Scenario #2

sc = 2; %scenario

aNA = 0.4;

%proportion of disputes between N & A when A is chosen

nNA = 1- aNA;

%proportion of disputes between N & M when N is chosen

nNM = 0.8;

%propotion of disputes between N & M when M is chosen

mNM = 1 - nNM;

%propotion of disputes between M & A when M is chosen

mMA = 0.3;

%propotion of disputes between M & A when A is chosen

aMA = 1 - mMA;

%proportion of population that has disputes when randomly mixed

k = 0.05;

%time step

delt = 1; %months

T = 300; %number of months to simulate

%% Scenario #3

sc = 3;

aNA = 0.75;

%proportion of disputes between N & A when A is chosen

nNA = 1-aNA;

%proportion of disputes between N & M when N is chosen

nNM = 0.4;

%propotion of disputes between N & M when M is chosen

mNM = 1 - nNM;

%propotion of disputes between M & A when M is chosen

mMA = 0.4;

%propotion of disputes between M & A when A is chosen

aMA = 1 - mMA;

%proportion of population that has disputes when randomly mixed

k = 0.3;

%time step

delt = 1; %months

T = 100; %number of months to simulate

%% Initial conditions

clear N A M

%initial proportion of negotiation

N(1,1) = 500;

%arbitration

A(1,1) = 500;

%mediation

M(1,1) = 500;

%total population

Ptot = N(1)+A(1)+M(1);

%% solve system over time using Euler's method

for j = 1:T

N(j+1,1) = (k/Ptot)*((nNA-aNA)*N(j,1)*A(j,1) + (nNM-mNM)*N(j,1)*M(j,1))+N(j,1);

A(j+1,1) = (k/Ptot)*((aNA-nNA)*N(j,1)*A(j,1) + (aMA-mMA)*A(j,1)*M(j,1))+A(j,1);

M(j+1,1) = (k/Ptot)*((mNM-nNM)*N(j,1)*M(j,1) + (mMA-aMA)*A(j,1)*M(j,1))+M(j,1);

end

subplot(1,3,sc)

hold all

plot((1:T+1),N,'linewidth',2)

plot((1:T+1),A,'linewidth',2)

plot((1:T+1),M,'linewidth',2)

legend('Negotiators','Arbitrators','Mediators')

title(['Scenario #',num2str(sc)],'fontsize',14)

xlabel('Time (months)')

xlim([0 T])

ylabel('Population')

set(gca,'fontsize',14)

grid on

box on

%check sum of all three is constant

check = N+A+M;

%% Truncation Error from Euler's

%initiate to 0

teN(1) = 0; %negotiator

teM(1) = 0; %mediator

teA(1) = 0; %arbitrator

for j = 1:T

%simulate 1st

N(j+1,1) = (k/Ptot)*((nNA-aNA)*N(j,1)*A(j,1) + (nNM-mNM)*N(j,1)*M(j,1))+N(j,1);

A(j+1,1) = (k/Ptot)*((aNA-nNA)*N(j,1)*A(j,1) + (aMA-mMA)*A(j,1)*M(j,1))+A(j,1);

M(j+1,1) = (k/Ptot)*((mNM-nNM)*N(j,1)*M(j,1) + (mMA-aMA)*A(j,1)*M(j,1))+M(j,1);

%calculate 2nd deriv at j+1

d2N(j+1,1) = (k/Ptot)*((nNA-aNA)*A(j,1) + (nNM-mNM)*M(j,1));

d2A(j+1,1) = (k/Ptot)*((aNA-nNA)*N(j,1) + (aMA-mMA)*M(j,1));

d2M(j+1,1) = (k/Ptot)*((mNM-nNM)*M(j,1) + (mMA-aMA)*M(j,1));

teN(j+1,1) = 0.5*d2N(j+1)*1^2;

teA(j+1,1) = 0.5*d2A(j+1)*1^2;

teM(j+1,1) = 0.5*d2M(j+1)*1^2;

end

%might need to add next order of taylor series

%plot

subplot(1,3,sc)

hold all

plot((1:T+1),teN,'linewidth',2)

plot((1:T+1),teA,'linewidth',2)

plot((1:T+1),teM,'linewidth',2)

legend('Negotiators','Arbitrators','Mediators')

title('Truncation Error','fontsize',14)

xlabel('Time (months)')

ylabel('Population')

set(gca,'fontsize',14)

grid on

box on

%% compare to alternative method - 5th order RK

clear Ark Nrk Mrk

%initial proportion of negotiation

Nrk(1,1) = 500;

%arbitration

Ark(1,1) = 500;

%mediation

Mrk(1,1) = 500;

for j = 1:T

%calculate k1 for N, A & M

k1N = fN(Nrk(j,1),Ark(j,1),Mrk(j,1),k,aNA,nNM,mMA);

k1A = fA(Nrk(j,1),Ark(j,1),Mrk(j,1),k,aNA,nNM,mMA);

k1M = fM(Nrk(j,1),Ark(j,1),Mrk(j,1),k,aNA,nNM,mMA);

%k2

k2N = fN(Nrk(j,1)+0.25*k1N,Ark(j,1)+(0.25)*k1A,...

Mrk(j,1)+0.25*k1M,k,aNA,nNM,mMA);

k2A = fA(Nrk(j,1)+0.25*k1N,Ark(j,1)+0.25*k1A,...

Mrk(j,1)+0.25*k1M,k,aNA,nNM,mMA);

k2M = fM(Nrk(j,1)+0.25*k1N,Ark(j,1)+0.25*k1A,...

Mrk(j,1)+0.25*k1M,k,aNA,nNM,mMA);

%k3

k3N = fN(Nrk(j,1)+(1/8)*k1N+(1/8)*k2N,Ark(j,1)+(1/8)*k1A+(1/8)*k2A,...

Mrk(j,1)+(1/8)*k1M+(1/8)*k2M,k,aNA,nNM,mMA);

k3A = fA(Nrk(j,1)+(1/8)*k1N+(1/8)*k2N,Ark(j,1)+(1/8)*k1A+(1/8)*k2A,...

Mrk(j,1)+(1/8)*k1M+(1/8)*k2M,k,aNA,nNM,mMA);

k3M = fM(Nrk(j,1)+(1/8)*k1N+(1/8)*k2N,Ark(j,1)+(1/8)*k1A+(1/8)*k2A,...

Mrk(j,1)+(1/8)*k1M+(1/8)*k2M,k,aNA,nNM,mMA);

%k4

k4N = fN(Nrk(j,1)-0.5*k2N+k3N,Ark(j,1)-0.5*k2A+k3A,...

Mrk(j,1)-0.5*k2M+k3M,k,aNA,nNM,mMA);

k4A = fA(Nrk(j,1)-0.5*k2N+k3N,Ark(j,1)-0.5*k2A+k3A,...

Mrk(j,1)-0.5*k2M+k3M,k,aNA,nNM,mMA);

k4M = fM(Nrk(j,1)-0.5*k2N+k3N,Ark(j,1)-0.5*k2A+k3A,...

Mrk(j,1)-0.5*k2M+k3M,k,aNA,nNM,mMA);

%k5

k5N = fN(Nrk(j,1)+(3/16)*k1N+(9/16)*k4N,...

Ark(j,1)+(3/16)*k1A+(9/16)*k4A,...

Mrk(j,1)+(3/16)*k1M+(9/16)*k4M,k,aNA,nNM,mMA);

k5A = fA(Nrk(j,1)+(3/16)*k1N+(9/16)*k4N,...

Ark(j,1)+(3/16)*k1A+(9/16)*k4A,...

Mrk(j,1)+(3/16)*k1M+(9/16)*k4M,k,aNA,nNM,mMA);

k5M = fM(Nrk(j,1)+(3/16)*k1N+(9/16)*k4N,...

Ark(j,1)+(3/16)*k1A+(9/16)*k4A,...

Mrk(j,1)+(3/16)*k1M+(9/16)*k4M,k,aNA,nNM,mMA);

%k6

k6N = fN(Nrk(j,1)-(3/7)*k1N+(2/7)*k2N+(12/7)*k3N-(12/7)*k4N+(8/7)*k5N,...

Ark(j,1)-(3/7)*k1A+(2/7)*k2A+(12/7)*k3A-(12/7)*k4A+(8/7)*k5A,...

Mrk(j,1)-(3/7)*k1M+(2/7)*k2M+(12/7)*k3M-(12/7)*k4M+(8/7)*k5M,k,aNA,nNM,mMA);

k6A = fA(Nrk(j,1)-(3/7)*k1N+(2/7)*k2N+(12/7)*k3N-(12/7)*k4N+(8/7)*k5N,...

Ark(j,1)-(3/7)*k1A+(2/7)*k2A+(12/7)*k3A-(12/7)*k4A+(8/7)*k5A,...

Mrk(j,1)-(3/7)*k1M+(2/7)*k2M+(12/7)*k3M-(12/7)*k4M+(8/7)*k5M,k,aNA,nNM,mMA);

k6M = fM(Nrk(j,1)-(3/7)*k1N+(2/7)*k2N+(12/7)*k3N-(12/7)*k4N+(8/7)*k5N,...

Ark(j,1)-(3/7)*k1A+(2/7)*k2A+(12/7)*k3A-(12/7)*k4A+(8/7)*k5A,...

Mrk(j,1)-(3/7)*k1M+(2/7)*k2M+(12/7)*k3M-(12/7)*k4M+(8/7)*k5M,k,aNA,nNM,mMA);

%sumulate N

Nrk(j+1,1) = Nrk(j,1)+(1/90)*(7*k1N+32*k3N+12*k4N+32*k5N+7*k6N);

%simulate A

Ark(j+1,1) = Ark(j,1)+(1/90)*(7*k1A+32*k3A+12*k4A+32*k5A+7*k6A);

%sumulate M

Mrk(j+1,1) = Mrk(j,1)+(1/90)*(7*k1M+32*k3M+12*k4M+32*k5M+7*k6M);

end

%% plot results

figure

subplot(1,3,sc)

hold all

plot((1:T+1),Nrk,'linewidth',2)

plot((1:T+1),Ark,'linewidth',2)

plot((1:T+1),Mrk,'linewidth',2)

legend('Negotiators','Arbitrators','Mediators')

title(['Scenario #',num2str(sc)],'fontsize',14)

xlabel('Time (months)')

ylabel('Population')

set(gca,'fontsize',14)

grid on

box on

xlim([0 T+1])

%check sum = 1500

check2 = Nrk+Ark+Mrk;

%% plot Euler's and Butchers together

figure

hold all

plot((1:T+1),N,'linewidth',2,'color',[0,0.447,0.741])

plot((1:T+1),Nrk,'--','linewidth',2,'color',[0,0.447,0.741])

plot((1:T+1),A,'linewidth',2,'color',[0.85,0.325,0.098])

plot((1:T+1),Ark,'--','linewidth',2,'color',[0.85,0.325,0.098])

plot((1:T+1),M,'linewidth',2,'color',[0.929,0.694,0.125])

plot((1:T+1),Mrk,'--','linewidth',2,'color',[0.929,0.694,0.125])

legend('Neg - Euler''s','Neg - Butcher''s',...

'Arb - Euler''s','Arb - Butcher''s',...

'Med - Euler''s','Med - Butcher''s')

title('Comparison between Euler''s & Butcher''s Methods','fontsize',14)

xlabel('Time (months)')

ylabel('Differences in Populations')

set(gca,'fontsize',14)

grid on

box on

xlim([1 100])

%% Plot truncation error

figure

hold all

plot((1:T+1),abs(teN),'linewidth',2,'color',[0,0.447,0.741])

plot((1:T+1),abs(teA),'linewidth',2,'color',[0.85,0.325,0.098])

plot((1:T+1),abs(teM),'linewidth',2,'color',[0.929,0.694,0.125])

legend('Negotiators','Arbitrators',...

'Mediators')

title('Truncation Error from Euler''s Method','fontsize',14)

xlabel('Time (months)')

ylabel('Differences in Populations')

set(gca,'fontsize',14)

grid on

box on

xlim([1 100])

%% Plot difference between the two methods and compare to truncation error

%difference in numerical methods

Ndiff = abs(Nrk-N);

Adiff = abs(Ark-A);

Mdiff = abs(Mrk-M);

%plot differences with truncation error

figure

hold all

plot((1:T+1),abs(teN),'linewidth',2,'color',[0,0.447,0.741])

plot((1:T+1),Ndiff,'--','linewidth',2,'color',[0,0.447,0.741])

plot((1:T+1),abs(teA),'linewidth',2,'color',[0.85,0.325,0.098])

plot((1:T+1),Adiff,'--','linewidth',2,'color',[0.85,0.325,0.098])

plot((1:T+1),abs(teM),'linewidth',2,'color',[0.929,0.694,0.125])

plot((1:T+1),Mdiff,'--','linewidth',2,'color',[0.929,0.694,0.125])

legend('Neg - Truncation','Neg - Diff',...

'Arb - Truncation','Arb - Diff',...

'Med - Truncation','Med - Diff')

title('Comparison of Truncation Error to Difference in Numerical Methods','fontsize',14)

xlabel('Time (months)')

ylabel('Differences in Populations')

set(gca,'fontsize',14)

grid on

box on

xlim([1 100])

%% Fit population growth data

%load data

load('PopulationGrowth.mat')

%scatter plots of data

figure

%Negotiators

subplot(1,3,1)

scatter(num_month,Neg_obs,20,'*','markeredgecolor',[0,0.447,0.741])

title('Negotiators','fontsize',14)

xlabel('Time (months)')

ylabel('Incoming Populations')

set(gca,'fontsize',14)

grid on

box on

xlim([1 100])

%Mediators

subplot(1,3,2)

scatter(num_month,Med_obs,20,'*','markeredgecolor',[0.929,0.694,0.125])

title('Mediators','fontsize',14)

xlabel('Time (months)')

ylabel('Incoming Populations')

set(gca,'fontsize',14)

grid on

box on

xlim([1 100])

%Arbitrators

subplot(1,3,3)

scatter(num_month,Arb_obs,20,'*','markeredgecolor',[0.85,0.325,0.098])

title('Arbitrators','fontsize',14)

xlabel('Time (months)')

ylabel('Incoming Populations')

set(gca,'fontsize',14)

grid on

box on

xlim([1 100])

%% Linear model for negotiators

%build Z matrix

Z = ones(100,2);

Z(:,2) = (1:100); %number of months that have passed

%Y vector

Y = Neg_obs(1:100);

%solve for coefficients

A = inv(Z'*Z)*(Z'*Y);

a0 = A(1);

a1 = A(2);

%fitted mode

NegLin = a0+a1.*Z(:,2);

%model performance metrics

%LINEAR

% total sum of squares

St = sum((Neg_obs - mean(Neg_obs)).^2);

% sum of squared residuals

Sr = sum((Neg_obs - NegLin).^2);

% R-squared

r2 = (St - Sr)./St; %0.8563

%standard error of y given x

n = 100;

Syx = sqrt(Sr/(n-2)); %3.5181

%plot data versus model

figure

hold all

scatter(num_month,Neg_obs,20,'*','markeredgecolor',[0,0.447,0.741])

plot(num_month,NegLin,'-k','linewidth',1.5)

%add error

plot(num_month,NegLin+Syx,'--k','linewidth',1.5)

plot(num_month,NegLin-Syx,'--k','linewidth',1.5)

legend('Data','Linear Model','Standard Error')

grid on

title('Negotiators - Linear Model')

ylabel('Population Growth')

set(gca,'fontsize',14)

grid on

xlim([1 100])

box on

%% Saturation growth model for mediators

close all

%y=a3(x/()b3+x))

%linear form

%1/y=1/a3+(b3/a3)(1/x)

%build Z matrix

Z = ones(100,2);

Z(:,2) = 1./(1:100)'; %inverse of num_months

%Y vector

Y = 1./Med_obs;

%plot 1/y versus 1/x

figure

scatter(Y,Z(:,2),50,'*','markeredgecolor',[0.929,0.694,0.125])

title('Scatter Plot of Transformed Data')

ylabel('1/y')

xlabel('1/x')

set(gca,'fontsize',14)

grid on

box on

%I think 1st value is throwing this off

%ignore 1st value

Z(1,:) = [];

Y(1) = [];

%solve for coefficients

A = (Z'*Z)\(Z'*Y);

%solve for values in sat growth rate model

a3 = 1/A(1);

b3 = A(2)*a3;

%fitted model

MedSG = a3.*((1:100)'./(b3+(1:100)'));

%model performance metrics

% total sum of squares

St = sum((Med_obs - mean(Med_obs)).^2);

% sum of squared residuals

Sr = sum((Med_obs - MedSG).^2);

% R-squared

r2 = (St - Sr)./St;

%standard error of y given x

n = 100;

Syx = sqrt(Sr/(n-2));

%plot data versus model

figure

hold all

scatter(num_month,Med_obs,20,'*','markeredgecolor',[0.929,0.694,0.125])

plot(num_month,MedSG,'-k','linewidth',1.5)

%add error

plot(num_month,MedSG+Syx,'--k','linewidth',1.5)

plot(num_month,MedSG-Syx,'--k','linewidth',1.5)

legend('Data','Saturation Growth Model','Standard Error')

grid on

title('Mediators - Saturation Growth Model')

ylabel('Population Growth')

set(gca,'fontsize',14)

grid on

xlim([1 100])

box on

%% Sinusoid for Arbitrators

%model A(t) = a0+C1cos(wot+theta)

%linear transform A(t) = a0+ A1*cos(wot)+B1sin(wot)

%period - assume 36 months (3 years)

T = 36;

%angular frequency

wo = (2*pi)/T;

%build Z matrix

Z = ones(100,3);

Z(:,2) = cos(wo.*num_month);

Z(:,3) = sin(wo.*num_month);

%observations

Y = Arb_obs;

%solve for coefficients

A = (Z'*Z)\(Z'*Y);

%get model values

a0 = A(1);

c1 = sqrt(A(2)^2+A(3)^2);

theta = atan(-A(3)/A(2));

%fitted model

ArbSin = a0+c1.*cos(wo.*num_month-theta);

%model performance metrics

% total sum of squares

St = sum((Arb_obs - mean(Arb_obs)).^2);

% sum of squared residuals

Sr = sum((Arb_obs - ArbSin).^2);

% R-squared

r2 = (St - Sr)./St;

%standard error of y given x

n = 100;

Syx = sqrt(Sr/(n-2));

%plot data versus model

figure

hold all

scatter(num_month,Arb_obs,20,'*','markeredgecolor',[0.85,0.325,0.098])

plot(num_month,ArbSin,'-k','linewidth',1.5)

%add error

plot(num_month,ArbSin+Syx,'--k','linewidth',1.5)

plot(num_month,ArbSin-Syx,'--k','linewidth',1.5)

legend('Data','Sinusoid Model','Standard Error')

grid on

title('Arbitrators - Sinusoid Model')

ylabel('Population Growth')

set(gca,'fontsize',14)

grid on

xlim([1 100])

box on

%% Simulate model with population increase

T = 300;

tvec = (0:T)';

%Linear model for negotiators Nin = 5+0.3*t

Nin = 5.3469+0.2945.*tvec;

%saturation growth rate for mediators

Min = 28.387.*(tvec./(13.055+tvec));

%sinusoid for arbitrators

Ain = 14.8657+7.6007.*cos(0.1745.*tvec+1.5561);

figure

hold all

plot(tvec,Nin,'linewidth',2)

plot(tvec,Min,'linewidth',2)

plot(tvec,Ain,'linewidth',2)

legend('Incoming Negotiators','Incoming Mediators','Incoming Arbitrators')

xlabel('Time (months)')

ylabel('Incoming Population')

title('Population Increase Over Time')

set(gca,'fontsize',14)

grid on

box on

%have to calculate total population to get proportion of each group coming

%in

Ptot(1) = 1500; %initial population

for j = 1:T

Ptot(j+1,1) = Ptot(j,1)+Nin(j)+Min(j)+Ain(j);

end

%% Scenario #1

sc = 1; %scenario

%proportion of disputes between N & A when N is chosen as resoltuion

%strategy

aNA = 0.7;

%proportion of disputes between N & A when A is chosen

nNA = 1- aNA;

%proportion of disputes between N & M when N is chosen

nNM = 0.4;

%propotion of disputes between N & M when M is chosen

mNM = 1 - nNM;

%propotion of disputes between M & A when M is chosen

mMA = 0.8;

%propotion of disputes between M & A when A is chosen

aMA = 1 - mMA;

%proportion of population that has disputes when randomly mixed

k = 0.2;

%time step

delt = 1; %months

T = 100; %number of months to simulate

%% Scenario #2

sc = 2; %scenario

aNA = 0.4;

%proportion of disputes between N & A when A is chosen

nNA = 1- aNA;

%proportion of disputes between N & M when N is chosen

nNM = 0.8;

%propotion of disputes between N & M when M is chosen

mNM = 1 - nNM;

%propotion of disputes between M & A when M is chosen

mMA = 0.3;

%propotion of disputes between M & A when A is chosen

aMA = 1 - mMA;

%proportion of population that has disputes when randomly mixed

k = 0.05;

%time step

delt = 1; %months

T = 100; %number of months to simulate

%% Scenario #3

sc = 3;

aNA = 0.75;

%proportion of disputes between N & A when A is chosen

nNA = 1-aNA;

%proportion of disputes between N & M when N is chosen

nNM = 0.4;

%propotion of disputes between N & M when M is chosen

mNM = 1 - nNM;

%propotion of disputes between M & A when M is chosen

mMA = 0.4;

%propotion of disputes between M & A when A is chosen

aMA = 1 - mMA;

%proportion of population that has disputes when randomly mixed

k = 0.3;

%time step

delt = 1; %months

T = 100; %number of months to simulate

%% solve system again using Euler's with incoming population functions

clear N A M

%initial proportion of negotiation

N(1,1) = 500;

%arbitration

A(1,1) = 500;

%mediation

M(1,1) = 500;

for j = 1:T

N(j+1,1) = (k/Ptot(j,1))*((nNA-aNA)*N(j,1)*A(j,1) + (nNM-mNM)*N(j,1)*M(j,1))+N(j,1)...

+Nin(j,1);

A(j+1,1) = (k/Ptot(j,1))*((aNA-nNA)*N(j,1)*A(j,1) + (aMA-mMA)*A(j,1)*M(j,1))+A(j,1)...

+Ain(j,1);

M(j+1,1) = (k/Ptot(j,1))*((mNM-mNM)*N(j,1)*M(j,1) + (mMA-aMA)*A(j,1)*M(j,1))+M(j,1)...

+Min(j,1);

end

%% plot for each scenario

figure

%sc refers to scenario

subplot(1,3,sc)

hold all

plot((1:T+1),N,'linewidth',2)

plot((1:T+1),A,'linewidth',2)

plot((1:T+1),M,'linewidth',2)

legend('Negotiators','Arbitrators','Mediators')

title(['Scenario #',num2str(sc)],'fontsize',14)

xlabel('Time (months)')

ylabel('Population')

set(gca,'fontsize',14)

grid on

box on

%% Simulate model with regression errors

%use standard error to do 68% confidence interval

Nin1 = 5.3469+0.2945.*tvec+3.5181;

Nin2 = 5.3469+0.2945.*tvec-3.5181;

%saturation growth rate for mediators

Min1 = 28.387.*(tvec./(13.055+tvec))+0.9937;

Min2 = 28.387.*(tvec./(13.055+tvec))-0.9937;

%sinusoid for arbitrators

Ain1 = 14.8657+7.6007.*cos(0.1745.*tvec+1.5561)+2.4618;

Ain2 = 14.8657+7.6007.*cos(0.1745.*tvec+1.5561)-2.4618;

clear N A M

%initial proportion of negotiation

N1(1,1) = 500; N2(1,1) = 500;

%arbitration

A1(1,1) = 500; A2(1,1) = 500;

%mediation

M1(1,1) = 500; M2(1,1) = 500;

for j = 1:T

N1(j+1,1) = (k/Ptot(j,1))*((nNA-aNA)*N1(j,1)*A1(j,1) + (nNM-mNM)*N1(j,1)*M1(j,1))+N1(j,1)...

+Nin1(j,1);

N2(j+1,1) = (k/Ptot(j,1))*((nNA-aNA)*N2(j,1)*A2(j,1) + (nNM-mNM)*N2(j,1)*M2(j,1))+N2(j,1)...

+Nin2(j,1);

A1(j+1,1) = (k/Ptot(j,1))*((aNA-nNA)*N1(j,1)*A1(j,1) + (aMA-mMA)*A1(j,1)*M1(j,1))+A1(j,1)...

+Ain1(j,1);

A2(j+1,1) = (k/Ptot(j,1))*((aNA-nNA)*N2(j,1)*A2(j,1) + (aMA-mMA)*A2(j,1)*M2(j,1))+A2(j,1)...

+Ain2(j,1);

M1(j+1,1) = (k/Ptot(j,1))*((mNM-mNM)*N1(j,1)*M1(j,1) + (mMA-aMA)*A1(j,1)*M1(j,1))+M1(j,1)...

+Min1(j,1);

M2(j+1,1) = (k/Ptot(j,1))*((mNM-mNM)*N2(j,1)*M2(j,1) + (mMA-aMA)*A2(j,1)*M2(j,1))+M2(j,1)...

+Min2(j,1);

end

%% plot for each scenario

figure

subplot(1,3,sc)

hold all

plot((1:T+1),N1,'--','linewidth',2,'Color',[0,0.447,0.741])

plot((1:T+1),A1,'--','linewidth',2,'color',[0.85,0.325,0.098])

plot((1:T+1),M1,'--','linewidth',2,'color',[0.929,0.694,0.125])

plot((1:T+1),N2,'--','linewidth',2,'Color',[0,0.447,0.741])

plot((1:T+1),A2,'--','linewidth',2,'color',[0.85,0.325,0.098])

plot((1:T+1),M2,'--','linewidth',2,'color',[0.929,0.694,0.125])

legend('Negotiators','Arbitrators','Mediators')

title(['Scenario #',num2str(sc)],'fontsize',14)

xlabel('Time (months)')

ylabel('Population')

set(gca,'fontsize',14)

grid on

box on

%% Operating costs

%cost per person per month

MonthCost = ([400 450 475 375 300 290 325 340 345 370 395 410])';

%make a longer vector

j=0;

for yr = 1:25

OperCost(j+1:j+12,1) = MonthCost;

j = length(OperCost);

end

%calculate total for each group - simplest form of numerical integration

Nopc = N.*OperCost(1:length(N));

Mopc = M.*OperCost(1:length(M));

Aopc = A.*OperCost(1:length(A));

%total operating costs

opcTot = Nopc+Mopc+Aopc;

figure

hold all

plot((1:T+1),Nopc,'linewidth',2)

plot((1:T+1),Aopc,'linewidth',2)

plot((1:T+1),Mopc,'linewidth',2)

plot((1:T+1),opcTot,'linewidth',2,'color',[0 0 0])

legend('Negotiators','Arbitrators','Mediators','Total')

title('Operating Costs over Time','fontsize',14)

xlabel('Time (months)')

ylabel('Total Cost ($)')

set(gca,'fontsize',14)

grid on

box on

xlim([1 100])

%% Trapezoid rule for operating costs

for j = 1:100

Ct(j) = 0.25*(Ptot(j)+Ptot(j+1))*(OperCost(j)+OperCost(j+1));

end

CtotTrap = sum(Ct);

%% Cost from damage and injuries

%monthly damage for neg, arb, med

DmgCost = [150 45 20];

InjCost = [120 25 50];

%easiest way to calculate (rectangles instead of trapezoids)

%calculate monthly costs from Damage for each group

Ndmg = N.*DmgCost(1);

Admg = A.*DmgCost(2);

Mdmg = M.*DmgCost(3);

%total Damage

dmgTot = Ndmg+Admg+Mdmg;

%calculate monthly costs from Injury for each group

Ninj = N.*InjCost(1);

Ainj = A.*InjCost(2);

Minj = M.*InjCost(3);

%total injury cost

injTot = Ninj+Ainj+Minj;

%% Apply trapezoid rule to compute monthly cost of damage and injury

for j = 1:100

Cd(j,1) = 0.5*(N(j)+N(j+1))*DmgCost(1)+0.5*(M(j)+M(j+1))*DmgCost(2)...

+0.5*(A(j)+A(j+1))*DmgCost(3);

Ci(j,1)= 0.5*(N(j)+N(j+1))*InjCost(1)+0.5*(M(j)+M(j+1))*InjCost(2)...

+0.5*(A(j)+A(j+1))*InjCost(3);

end

Ctotd = sum(Cd);

Ctoti = sum(Ci);

%% plot individual costs on separate plots

figure

%operating costs

subplot(3,1,1)

hold all

plot((1:T+1),Nopc,'linewidth',2)

plot((1:T+1),Aopc,'linewidth',2)

plot((1:T+1),Mopc,'linewidth',2)

legend('Negotiators','Arbitrators','Mediators')

title('Operating Costs over Time','fontsize',14)

xlabel('Time (months)')

ylabel('Total Cost ($)')

set(gca,'fontsize',14)

grid on

box on

xlim([1 100])

%damage costs

subplot(3,1,2)

hold all

plot((1:T+1),Ndmg,'linewidth',2)

plot((1:T+1),Admg,'linewidth',2)

plot((1:T+1),Mdmg,'linewidth',2)

legend('Negotiators','Arbitrators','Mediators')

title('Damge Costs over Time','fontsize',14)

xlabel('Time (months)')

ylabel('Total Cost ($)')

set(gca,'fontsize',14)

grid on

box on

xlim([1 100])

%injury costs

subplot(3,1,3)

hold all

plot((1:T+1),Ninj,'linewidth',2)

plot((1:T+1),Ainj,'linewidth',2)

plot((1:T+1),Minj,'linewidth',2)

legend('Negotiators','Arbitrators','Mediators')

title('Injury Costs over Time','fontsize',14)

xlabel('Time (months)')

ylabel('Total Cost ($)')

set(gca,'fontsize',14)

grid on

box on

xlim([1 100])

%% total from all three groups compared to one another

figure

hold all

plot((1:T+1),opcTot,'linewidth',2)

plot((1:T+1),dmgTot,'linewidth',2)

plot((1:T+1),injTot,'linewidth',2)

legend('Operating Costs','Damage Costs','Injury Costs')

title('Comparisons of Cost over Time','fontsize',14)

xlabel('Time (months)')

ylabel('Total Cost ($)')

set(gca,'fontsize',14)

grid on

box on

%% function for dN/dt

function out = fN(N,A,M,k,aNA,nNM,mMA)

nNA = 1- aNA;

mNM = 1 - nNM;

aMA = 1 - mMA;

Ptot = 1500;

out = (k/Ptot)*((nNA-aNA)*A*N + (nNM-mNM)*N*M);

end

%% function for dA/dt

function out = fA(N,A,M,k,aNA,nNM,mMA)

nNA = 1- aNA;

mNM = 1 - nNM;

aMA = 1 - mMA;

Ptot = 1500;

out = (k/Ptot)*((aNA-nNA)*A*N +(aMA-mMA)*A*M);

end

%% function for dM/dt

function out = fM(N,A,M,k,aNA,nNM,mMA)

nNA = 1- aNA;

mNM = 1 - nNM;

aMA = 1 - mMA;

Ptot = 1500;

out = (k/Ptot)*((mMA-aMA)*M*A +(mNM-nNM)*N*M);

end