JointMDPNonlinearEqn.run

# Load model and data
model "MarkovActv.mod";
data "MarkovActv.dat";

# Is Cauchy distribution used ? 
fix IS_CAUCHY := 1;

# logging options
let debug_log := 1;

# include the code that define the component state and choice set
include MDPStateAction.run
include JointMDPStateAction.run

# Define the composite MDP problem for the entire household
problem compositeMDP:

# Choose the objective function
likelihood0,

# List the variables
EW, 
actvUtil,
sumActvUtil,
jointChoiceUtil,
jointChoiceProb,
# lower, upper, 
# actvUtil, sumActvUtil, sumTravelCost, jointActvUtil, jointChoiceUtil, 
Um, b, c, rho,
# VoT, theta, 
# EV, actvUtil, choiceUtil, choiceProb, VoT, theta, Uw, xi, gamma,

# List the constraints
# jointChoiceProb_Range,
Bellman_Joint,
Bellman_JointH;
# Bellman_Lower,
# Bellman_LowerH,
# Bellman_Upper,
# Bellman_UpperH,
# LowerBound, 
# UpperBound;

fix {n in PERS, j in ACTV} Um[n,j] := Um0[n,j];
fix {n in PERS, j in ACTV}  b[n,j] := b0[n,j];
fix {n in PERS, j in ACTV}  c[n,j] := c0[n,j];
fix {j in ACTV} rho[j] := rho0[j];

# Initialize the lower and upper bounds and the EW values
# for {(t,j1,j2) in XX} {
# 	fix lower[t,j1,j2] := max(EV[1,t,j1], EV[2,t,j2]);
# 	fix upper[t,j1,j2] := EV[1,t,j1] + EV[2,t,j2];
# 	let EW[t,j1,j2] := .5 * (lower[t,j1,j2] + upper[t,j1,j2]);
# 	# fix lower[t,j1,j2] := -500;
# 	# fix upper[t,j1,j2] := 500;
# }


# Specify KNITRO solver options:
option solver "C:\Ziena\KNITRO900\knitroampl\knitroampl.exe";
option knitro_options "alg=2 hessopt=1 outlev=3 maxit=500 xtol=0.0000000001 wantsol=1";
# option presolve 0;

# Solve command
solve;

# Solution status
printf "Composite MDP\n" > DATA/jointMDP.sol;
if match (solve_message, "Locally optimal solution") > 0 then { 
	printf "%s\n", "Optimal solution found" > DATA/jointMDP.sol;
	printf "%5.0f \n", 0 > DATA/jointMDP.sol ; 
}
else if match (solve_message, "Iteration limit reached") > 0 then { 
	printf "%s\n", "Iteration limit reached" > DATA/jointMDP.sol;
	printf "%5.0f \n", 400 > DATA/jointMDP.sol ; 
}
else if match (solve_message, "Evaluation error") > 0 then { 
	printf "%s\n", "Evaluation error" > DATA/jointMDP.sol;
	printf "%5.0f \n", 502 > DATA/jointMDP.sol ; 
}
else {
	printf "%s\n", "No solution" > DATA/jointMDP.sol;
	printf "%5.0f \n", 1000 > DATA/jointMDP.sol ; 
}

# display _solve_time > DATA/jointMDP.sol;

# write the structural parameters
display beta, VoT, theta > DATA/jointMDP.sol;
display Um, b, c, rho > DATA/jointMDP.sol;

# also display the structural parameters
printf "\n\nComposite MDP";
printf "\n\nstructural parameters (all fixed):\n\n";
display beta, VoT, theta;
display Um, b, c, rho;
# display Uw, xi, gamma, lambda;


# Output commands
option display_round 6, display_width 120;

# write the activity utility and travel cost
# display actvUtil > DATA/actvUtil.dat;
# display sumActvUtil > DATA/sumActvUtil.dat;
# display sumTravelCost > DATA/sumTravelCost.dat;
# display jointActvUtil > DATA/jointActvUtil.dat;


# export flow variables
include ExportJointNetworkFlows.run;

# export EW to file jointEUtil.m
printf "EW = zeros(%d, %d, %d);\n", H, M, M > DATA/jointEUtil.m;
for {t in TIME} {
	printf "EW(%3d, :, :) = [", t+1 > DATA/jointEUtil.m;
	for {j1 in AUW[1]} {
		for {j2 in AUW[2]} {
			if (t,j1,j2) in XX then
				printf " %8.2f", EW[t,j1,j2] > DATA/jointEUtil.m;
			else
				printf " %8s", "nan" > DATA/jointEUtil.m;
		}
		printf ";" > DATA/jointEUtil.m;
	}
	printf "];\n" > DATA/jointEUtil.m;
}

# start debug log
if debug_log == 1 then {

# export choiceProb to file jointChoiceProb.m
printf "Pj = zeros(%d, %d, %d, %d, %d, %d);\n", H, M, M, M, M, DH > DATA/jointChoiceProb.m;
for {(t,j1,j2) in XX} {
	for {h in 1..DH} {
		printf "Pj(%3d, %d, %d, :, :, %2d) = [", t+1, j1, j2, h > DATA/jointChoiceProb.m;
		for {a1 in AUW[1]} {
			for {a2 in AUW[2]} {
				if (a1,a2,h) in DD[t,j1,j2] then
					if jointChoiceProb[t,j1,j2,a1,a2,h] >= -0.000001 and 
					   jointChoiceProb[t,j1,j2,a1,a2,h] <=  1.000001 then
						printf " %8.5f", jointChoiceProb[t,j1,j2,a1,a2,h] > DATA/jointChoiceProb.m;
					else
						printf " %8s", "nan" > DATA/jointChoiceProb.m;
				else
					printf " %8.5f", 0.0 > DATA/jointChoiceProb.m;
			}
			printf ";" > DATA/jointChoiceProb.m;
		}
		printf "];\n" > DATA/jointChoiceProb.m;
	}
}

# export exp(jointChoiceProb) to file expJointChoiceProb.txt
printf "ePj = zeros(%d, %d, %d, %d, %d, %d);\n", H, M, M, M, M, DH > DATA/expJointChoiceProb.txt;
for {(t,j1,j2) in XX} {
	for {h in 1..DH} {
		printf "ePj(%3d, %d, %d, :, :, %2d) = [", t+1, j1, j2, h > DATA/expJointChoiceProb.txt;
		for {a1 in AUW[1]} {
			for {a2 in AUW[2]} {
				if (a1,a2,h) in DD[t,j1,j2] then
					printf " %8.2f", ( theta * ( jointChoiceUtil[t,j1,j2,a1,a2,h] + 
				  								 beta**h * EW[t,a1,a2]) - 
		 							   theta * EW[t,j1,j2]) > DATA/expJointChoiceProb.txt;
				else
					printf " %8s", "nan" > DATA/expJointChoiceProb.txt;
			}
			printf ";" > DATA/expJointChoiceProb.txt;
		}
		printf "];\n" > DATA/expJointChoiceProb.txt;
	}
}

# export jointChoiceUtil + beta^h*EW to file jointChoiceUtil&EW.txt
printf "ePj = zeros(%d, %d, %d, %d, %d, %d);\n", H, M, M, M, M, DH > DATA/jointChoiceUtil&EW.txt;
for {(t,j1,j2) in XX} {
	for {h in 1..DH} {
		printf "ePj(%3d, %d, %d, :, :, %2d) = [", t+1, j1, j2, h > DATA/jointChoiceUtil&EW.txt;
		for {a1 in AUW[1]} {
			for {a2 in AUW[2]} {
				if (a1,a2,h) in DD[t,j1,j2] then
					printf " %8.2f + %8.2f", jointChoiceUtil[t,j1,j2,a1,a2,h],
				  							 beta**h * EW[t,a1,a2] > DATA/jointChoiceUtil&EW.txt;
				else
					printf " %8s + %8.2f", "nan", beta**h * EW[t,a1,a2] > DATA/jointChoiceUtil&EW.txt;
			}
			printf ";" > DATA/jointChoiceUtil&EW.txt;
		}
		printf "];\n" > DATA/jointChoiceUtil&EW.txt;
	}
}

# end of log
}

# export fixed parameters to jointMDP.m
printf "beta = %f;\n", beta > DATA/jointMDP.m;
printf "theta = %f;\n", theta > DATA/jointMDP.m;

# export the structural parameters to jointMDP.m
printf "VoT0 = %f;\n", VoT > DATA/jointMDP.m;
printf "Um0 = [ " > DATA/jointMDP.m;
for {n in PERS} {
	for {j in ALLACTV} {
		printf "%f ", Um[n,j] > DATA/jointMDP.m;
	}
	printf "\n" > DATA/jointMDP.m;
}
printf "];\n" > DATA/jointMDP.m;
printf "b0 = [ " > DATA/jointMDP.m;
for {n in PERS} {
	for {j in ALLACTV} {
		printf "%f ", b[n,j] > DATA/jointMDP.m;
	}
	printf "\n" > DATA/jointMDP.m;
}
printf "];\n" > DATA/jointMDP.m;
printf "c0 = [ " > DATA/jointMDP.m;
for {n in PERS} {
	for {j in ALLACTV} {
		printf "%f ", c[n,j] > DATA/jointMDP.m;
	}
	printf "\n" > DATA/jointMDP.m;
}
printf "];\n" > DATA/jointMDP.m;