# ICE Project: "An Optimal Rule of Thumb for Pollution Permits Allocation"  
# by Mar Reguant and Eva Dardati
# August, 2008
#
# Questions or comments to mreguant@mit.edu
#
# References: Please, have a look at Sven Leyffer, Todd Munson and Karl Schmedders tutorials
# Model based on "Cournot Game with Learning and Investment" in K. Schmedders slides and
# "Incorporating oligopoly, CO2 emissions trading and green certificates into a
# power generation expansion model" by Pedro Linares et al, Automatica.

reset;
#option log_file "permit_path.log";
# option log_file "permit_knitro.log";

# choose the solver
# PATH
model permit_path.mod;
option presolve 0;
option solver pathampl;
solve;

# KNITRO
# model permit_knitro.mod;
# option presolve 0;
# option solver knitroampl;
# solve;

printf "" > results.out;
for{i in firms, (v,w) in states} {
	printf "%d , %d , %d , %f , %f , %f , \n", i, v, w, q[i,v,w], x[i,v,w], Pi[i,v,w] >> results.out;
}

# DISPLAY
# display q, x, z;
# display Pr;
# display V;

# COMPUTE LONG RUN DISTRIBUTIONS
# param init := nStates;
# param T := 100000;
# param theta_dis{firms cross 1..T}, integer;
# param q_dis{1..T};
# param e_dis{1..T};

# for the loop
# param thetaold{firms}, default init;
# param thetanew{firms};
# param u{firms cross 1..T} := Uniform(0,1);

# for {t in 1 .. T} {

	# let{i in firms} thetanew[i] := 
		# (if u[i,t] <= Pr[i, thetaold[i], thetaold[nextw(i)], -1] then thetaold[i] - 1) +
		# (if u[i,t] > Pr[i, thetaold[i], thetaold[nextw(i)], -1] and
		# u[i,t] <= Pr[i, thetaold[i], thetaold[nextw(i)], -1] + Pr[i, thetaold[i], thetaold[nextw(i)], 0] then thetaold[i]) +
		# (if  u[i,t] > 1-Pr[i, thetaold[i], thetaold[nextw(i)], 1] then thetaold[i] + 1);
	
	# let q_dis[t] := sum{i in firms} q[i,thetaold[i], thetaold[nextw(i)]];
	# let e_dis[t] := sum{i in firms} theta[i,thetaold[i], thetaold[nextw(i)]]*q[i,thetaold[i],thetaold[nextw(i)]];
	
	# let{i in firms} theta_dis[i,t] := thetanew[i];
	# let{i in firms} thetaold[i] := thetanew[i];
# }

# param theta_freq{i in firms, s in 1..nStates} :=		# distribution of theta
	# sum{t in 1..T} (if theta_dis[i,t] = s then 1)/T;
# printf "" > dist.out;
# for{i in firms, s in 1 .. nStates} {
	# printf "%d , %d , %f , \n", i, s, theta_freq[i,s] >> dist.out;
# }

# printf "" > qdist.out;						# vector of quantities and emissions
# printf "" > edist.out;
# for{t in 1..T} {
	# printf "%d , %f \n", t, q_dis[t] >> qdist.out;
	# printf "%d , %f \n", t, e_dis[t] >> edist.out;
# }