***********************************************************************************;
* RUN OPTIMIZATION MODEL
***********************************************************************************;

proc optmodel &presolver_options.;

	*Define data elements;
	put 'optmodel: Define data elements';
	set <num> ORDERS;
	set <num> LOCATIONS;
	num distance {ORDERS, LOCATIONS} init &inf_dist.;

	*Define variable 1;
	var on {LOCATIONS} binary init 0;

	*Read datasets 1;
	put 'optmodel: Read datasets 1';
	read data test.orders    into ORDERS=[ORDER_ID];
	read data test.points    into LOCATIONS=[LOCATION_ID] on=on;
	read data test.distances into [ORDER_ID LOCATION_ID] distance=&distance_var.;

	*Define variable 2;
	set VALID_MAPS = {p in ORDERS, w in LOCATIONS: distance[p,w] <= (&inf_dist.-10)};
	var map {VALID_MAPS} binary init 0;

	*Read datasets 2;
	put 'optmodel: Read datasets 2';
	read data test.distances into [ORDER_ID LOCATION_ID] map=map;

	*Set priority and branching direction (1=Round up to nearest integer);
	for {w in LOCATIONS}      on[w].priority = 5;
	for {w in LOCATIONS}      on[w].direction = 1;
	for {<p,w> in VALID_MAPS} map[p,w].direction = 1;

	*Calculations;
	put 'optmodel: Calculations';
	impvar total_LOCATIONs_on = sum{w in LOCATIONS} on[w];
	impvar total_maps_on  = sum{<p,w> in VALID_MAPS} map[p,w];
	impvar total_distance = sum{<p,w> in VALID_MAPS} map[p,w]*distance[p,w];
	impvar total_penalty  = total_LOCATIONs_on * &LOCATION_distance_penalty.;

	*Mechanical constraints;
	put 'optmodel: Mechanical constraints';
	con ORDERs_mapped_correctly    {p in ORDERS}:    (sum{<(p),w> in VALID_MAPS} map[p,w]) =  1;
	con LOCATIONs_mapped_correctly {w in LOCATIONS}: (sum{<p,(w)> in VALID_MAPS} map[p,w]) <= on[w] * &max_ORDERs_per_LOCATION.;

	*Business constraints;
	put 'optmodel: Business constraints';
	con max_LOCATIONs_on:  total_LOCATIONs_on  <= &final_max_LOCATIONs_on.;

	*Minimize distance;
	put 'optmodel: Set Goals';
	min goal = total_distance + total_penalty;
	save mps test.model_result_mps;

	*Initial solution check (seed);
	print "VALUES BEFORE SOLVING" total_LOCATIONs_on total_maps_on goal total_distance total_penalty;

	*Solve problem;
	put 'optmodel: Solve problem';
	/*solve with MILP / &milp_config.; */
	save mps MPSDATA;
	create data PRIMALIN(drop=j) from [j] _VAR_=_VAR_.label _VALUE_=_VAR_.sol;
	submit;
	proc optmilp 
		data=MPSDATA
		primalin=PRIMALIN
		primalout=PRIMALOUT 
		&milp_config.;
	run;
	endsubmit;
	read data PRIMALOUT into [_N_] _VAR_[_N_]=_VALUE_;

	*Print results;
	print "VALUES AFTER SOLVING" total_LOCATIONs_on total_maps_on goal total_distance total_penalty;

	*Save results;
	put 'optmodel: Save results';
	create data test.model_result_legs      (where=(map>0.01)) from [ORDER_ID LOCATION_ID] = {ORDERS, LOCATIONS} map ;
	create data test.model_result_LOCATIONs (where=(on>0.01))  from [LOCATION_ID]          = {LOCATIONS} on;
	create data test.model_result_ORDERs                       from [ORDER_ID]             = {ORDERS};
run;