%display nonZeros
%data

/* PROBLEM

Simplified model of grain transport to develop proof of concept.

Producing regions vary acreage of crop under production to produce tonnage (seasonal, by grain). Grain yield per acre is constant across seasons. Producing regions have associated export ports and storage facilities. In a given season, for each producing region, grain is moved to either the export port or storage. Consuming regions have a single associated import port, and a demand by grain (not seasonal, i.e. constant across seasons). Ocean shipping routes are defined between export and import ports.  At present, there is only a single route between each port pair

Cost elements are:
	ocean shipping cost
	storage cost
	production cost

The objective is to minimize total cost over the seasons of the year

Proposed features/limits are:
	1 grain
	3 producing regions
	3 storage facilities (one per production region)
	4 seasons
	7 consuming regions
	1 export port per producing region
	1 import port per consuming region
	1 ocean transfer route between producing and consuming regions
	1 storage element per producing region
	2 route types (canal and non-canal)
	21 possible routes (some canal and some non-canal)
	Canal tolling depending on route type
	No river systems/river reaches
	No direct transfers
	No transport capacity constraints for port, reach, route, transfer to/from storage
	Capacity constraints on production acreage, storage

Producing regions
BrazilSouth & Argentina (South America), US Gulf Coast (Midwest), US Pacific Northwest

Consuming Regions
China, Europe, Mexico, Japan, Southeast Asia, FSU-MiddleEast, Korea, US, Brazil
*/

variables:

# PRODUCTION

# Area in production by region, season, and grain in hectares.
production_hectares[production_regions,seasons,grains]:real[0..];

# Production yield in metric tons by region, season, and grain.
production_mtons[production_regions,seasons,grains]:real[0..];

# Production cost in USD by region, season, and grain
production_cost[production_regions,seasons,grains]:real[0..];

# Sum of production costs in USD across regions, seasons, and grains
total_production_cost:real[0..];

# STORAGE

# Units added (positive) or removed (negative) from storage at the end of the season by region and season in metric tons.
seasonal_net_storage_mtons[production_regions,seasons,grains]:real;

# Units in storage throughout the whole season by region and season in metric tons. (Excess goods produced in season 2 are added at the end of season 2/beginning of season 3.)
storage_mtons[production_regions,seasons,grains]:real[0..];

# Storage cost by region and season in USD. Storage units will be charged at the end of the season for goods that have spent a whole season in storage. Newly added storage is not charged at this time.
storage_cost[production_regions,seasons,grains]:real[0..];

# Total storage cost in USD.
total_storage_cost:real[0..];

# TRANSPORT

# Cost in USD by route, season, and grain for ocean transport.
transport_cost[export_ports,import_ports,route_types,seasons,grains]:real[0..]; 

# Sum of transport cost in USD across route, season, and grain.
total_transport_cost:real[0..];

# Cost in USD by route, season, and grain for canal tolls.
toll_cost[export_ports,import_ports,route_types,seasons,grains]:real[0..]; 

# Sum of toll cost in USD across route, season, and grain.
total_toll_cost:real[0..];

# Route tonnage in metric tons by route, season, and grain.
route_mtons[export_ports,import_ports,route_types,seasons,grains]:real[0..];

# Export tonnage by port, season, and grain in metric tons.
export_mtons[export_ports,seasons,grains]:real[0..];

# Import tonnage by port, season, and grain in metric tons.
import_mtons[import_ports,seasons,grains]:real[0..];

# Determines intervals for linear approximation.
tonnage_intervals[seasons,grains]:real[0..]; #uji

# Consumption demand summed over whole world
worldwide_demand[seasons,grains]:real[0..];

# Points on the curve used for finding slope of interval
fx[export_ports,import_ports,seasons,grains,approximation_intervals]:real[0..];

# Slope of intervals
m[export_ports,import_ports,seasons,grains,approximation_intervals]:real[0..];

# Break x value down into where the value falls with respect to intervals
x_val_interval[export_ports,import_ports,route_types,seasons,grains,approximation_intervals]:real[0..]; #xji

# X value total for approximation
x_val[export_ports,import_ports,route_types,seasons,grains]:real[0..]; #xj

objectives:

# Minimize total cost
costs: total_production_cost + total_transport_cost + total_storage_cost + total_toll_cost -> min;

constraints:

# Add the x intervals together to get the total x (pg 574)
solve_x$1$ {
e in export_ports,
i in import_ports,
r in route_types,
s in seasons,
g in grains: 
x_val[e,i,r,s,g] = sum {
n in approximation_intervals: 
x_val_interval[e,i,r,s,g,n]
};
}
# Force x_(i+1) = 0 whenever x_i < u_i (pg 574)
special_condition$1${
e in export_ports,
i in import_ports,
r in route_types,
s in seasons,
g in grains,
n in approximation_intervals: 
{ n < len(approximation_intervals) && #Syntax error unexpected SYMBOL_VAR
x_val_interval[e,i,r,s,g,n] < tonnage_intervals[s,g,n]: 
x_val_interval[e,i,r,s,g,n+1] = 0;
}
}

# Sum demand to create upper limit for intervals
worldwide_demand$2$ {
s in seasons,
g in grains:
worldwide_demand[s,g] = sum {
c in consumption_regions:
consumption_demand_mtons[c,s,g]
};
}

# Create equal-sized intervals
tonnage_inter$1$ {
s in seasons,
g in grains:
tonnage_intervals[s,g] = worldwide_demand[s,g] * (1/len(approximation_intervals));
}

# Bound the x value by the possible x value for the interval
limit_x_1_$1$ {
e in export_ports,
i in import_ports,
r in route_types,
s in seasons,
g in grains,
n in approximation_intervals: 
x_val_interval[e,i,r,s,g,n] <= tonnage_intervals[s,g]; 
}

transport_cost$1$ {
e in export_ports,
i in import_ports,
r in route_types,
s in seasons,
g in grains: 
transport_cost[e,i,r,s,g] = sum{
n in approximation_intervals:
x_val_interval[e,i,r,s,g,n] * m[e,i,s,g,n] * route_nmiles[e,i,r] * transport_cost_per_mton_nmile #Syntax issue:  Product of variable has too much indexes
};
} 

#Get cost for the tonnage intervals
solve_fx$1$ {
e in export_ports,
i in import_ports,
s in seasons,
g in grains,
n in approximation_intervals:
fx[e,i,s,g,n] = exp(a[e,i]) * (n*tonnage_intervals[s,g])^b[e,i]; # Syntax issue with power
}

#Get slope for piece-wise approximation
solve_m$1$ {
e in export_ports,
i in import_ports,
s in seasons,
g in grains,
n in approximation_intervals: 
{
n = 1: 
m[e,i,s,g,n] = fx[e,i,s,g,n] * (1/tonnage_intervals[s,g]); # Arithmetic syntax issue
| n > 1: 
m[e,i,s,g,n] = (fx[e,i,s,g,n] - fx[e,i,s,g,n-1])/(n*tonnage_intervals[s,g] - (n-1)*tonnage_intervals[s,g,]); # Arithmetic syntax issue
}
}

# PRODUCTION

# Calculate production cost by region, season, and grain
prod_cost$1$ {
pr in production_regions,
s in seasons,
g in grains:
production_cost[pr,s,g] = production_hectares[pr,s,g] * production_cost_per_hectare[pr,g];
}

# Limit area in production for each region by the area available
production_hectares$1$ {
pr in production_regions,
s in seasons,
g in grains: 
production_hectares[pr,s,g] <= production_hectares_available[pr,s,g];
}

# Calculate total production cost
sum_of_production_cost: 
total_production_cost = sum{
pr in production_regions,
s in seasons, 
g in grains:
production_cost[pr,s,g]
};

# Calculate production tonnage by region, season, and grain
prod_tonnage$1$ {
pr in production_regions,
s in seasons,
g in grains: 
production_mtons[pr,s,g] = production_hectares[pr,s,g] * production_mtons_per_hectare[pr,g];
}

# STORAGE

# The amount stored needs to be limited by the full storage capacity.
limit_storage_units_capacity$1${
pr in production_regions, 
s in seasons,
g in grains:
storage_mtons[pr,s,g] <= storage_mtons_available[pr,g];
}

# Calculate the amount stored for each season.
calc_storage_units$1${
pr in production_regions, 
s in seasons,
g in grains:
{s = 1: 
storage_mtons[pr,s,g] = storage_mtons_initial[pr,g]; 
| s > 1: 
storage_mtons[pr,s,g] = storage_mtons[pr,s-1,g] + seasonal_net_storage_mtons[pr,s-1,g];
}
}

# Calculate the total cost of storage.
calc_total_storage_cost$1$: 
total_storage_cost = sum{
pr in production_regions,
s in seasons,
g in grains:
storage_cost[pr,s,g]
};

# Calculate the cost of storage.
calc_storage_cost$1${
pr in production_regions,
s in seasons,
g in grains:
storage_cost[pr,s,g] = storage_mtons[pr,s,g] * storage_cost_per_mton[pr,g];
}

# End with the same amount of storage that we start with.
set_end_storage{
pr in production_regions,
g in grains:
storage_mtons[pr,len(seasons),g] + seasonal_net_storage_mtons[pr,len(seasons),g] = storage_mtons_initial[pr,g];
}

# EXPORT

# Calculate tonnage exported by region, season, and grain
export_mtons$1$ {
pr in production_regions, 
s in seasons, 
g in grains: 
export_mtons[production_ports[pr],s,g] = production_mtons[pr,s,g] - seasonal_net_storage_mtons[pr,s,g];
}

# TRANSPORT

# Calculate cost of transportation
transport_cost$1$ {
e in export_ports,
i in import_ports,
s in seasons,
g in grains,
r in route_types: 
transport_cost[e,i,r,s,g] = route_mtons[e,i,r,s,g] * route_nmiles[e,i,r] * transport_cost_per_mton_nmile;
} 

# Calculate total cost of transportation
sum_of_transport_cost: 
total_transport_cost = sum{
e in export_ports,
i in import_ports,
s in seasons,
g in grains,
r in route_types:
transport_cost[e,i,r,s,g]
};

# Calculate cost of toll
toll_cost$1$ {
e in export_ports,
i in import_ports,
s in seasons,
g in grains,
r in route_types: 
toll_cost[e,i,r,s,g] = route_mtons[e,i,r,s,g] * toll_per_mton[r];
} 

# Calculate total cost of toll
sum_of_toll_cost: 
total_toll_cost = sum{
e in export_ports,
i in import_ports,
s in seasons,
g in grains,
r in route_types:
toll_cost[e,i,r,s,g]
};

# IMPORT

# Make sure import ports can satisfy consumption demand. No excess tonnage.
import_port_demand$1$ {
i in import_ports,
s in seasons, 
g in grains: 
import_mtons[i,s,g] = consumption_demand_mtons[import_regions[i],s,g];
}

# Make sure that import tonnage is met by sum of route tonnage.
imp_tonnage$2$ {
i in import_ports,
s in seasons, 
g in grains:  
import_mtons[i,s,g] = sum { 
e in export_ports,
r in route_types: 
route_mtons[e,i,r,s,g]
};
}

# Calculate export tonnage based on route tonnage
exp_tonnage$1$ {
e in export_ports,
s in seasons, 
g in grains:  
export_mtons[e,s,g] = sum { 
i in import_ports,
r in route_types: 
route_mtons[e,i,r,s,g]
};
}

# Calculate production tonnage based on route tonnage and storage tonnage
prod_tonnage$2$ {
pr in production_regions,
s in seasons, 
g in grains:  
production_mtons[pr,s,g] = sum { 
i in import_ports,
r in route_types:
route_mtons[production_ports[pr],i,r,s,g]
} + seasonal_net_storage_mtons[pr,s,g];
}