Solved: Re: Proc OptModel excessive solution times

genemroz · Posted 06-21-2021 11:36 AM

This is a follow-up question relate to my posting of yesterday. Rob Pratt responded with an excellent solution to my question that works well within the "basic" version that was posted with the question. But when I tried to implement his solution into a model version with more cameras, orientations and larger target areas, solution time for PROC OptModel becomes excessive. I didn't experience that problem with my earlier version of the optimization model. What's the best way to proceed?

Thanks,

Gene

RobPratt · Posted 06-21-2021 05:13 PM

For Version B with 420, the MILP presolver is able to remove almost all of the variables and constraints. For Version B with 440, it is not able to do so, for some reason that I cannot yet explain. You can save the presolver some work in this case by replacing the explicit variable XIT and constraint XIT_CON1 with the following implicit variable:

	impvar XIT {target in TARGETS} =
      sum{station in stations, orientation in ORIENTATIONS} 
   		matrix[station, orientation, target] * Chi[station, orientation];

Remember to declare the TriplyCoveredCon constraint after this IMPVAR statement.

With this change for Version B, both 420 and 440 (and even larger values of RoO) solve instantly.

View solution in original post

RobPratt · Posted 06-21-2021 01:01 PM

Please post the log and preferably also the data and code.

genemroz · Posted 06-21-2021 03:33 PM

Per Rob's request, I'm posting two versions of my camera network optimization code:

Version A (below) is the code (over 1000 lines) that includes my implementation of PROC OptModel prior to Rob Pratt's response to my request for a improved algorithm. The code is stand-alone and ready to run. At the top of the code, note the macro variable RoO. This macro variable sets the dimension of target area. It is set at 420. Run the code with RoO at 420 and then change it to 440 and run again. You'll see that the ProcOptModel solultion time is about 6 seconds in both cases. This is the control case.

Version B (also below) is the same code except that Proc OptModel has been modified to incorporate Rob Pratts' suggested code for the improvement I requested. This code is stand alone and ready to run also. As before, run the code with RoO at 420 and then change it to 440 and run again. You'll see that the Proc OptModel solution time increases from about 3 seconds to over 100 seconds when RoO is changed to 440. This seem "excessive" to me unless, perhaps, my expectations are unreasonable.

I would also note that the solution time gets much slower (over 5 minutes) as I increase RoO beyond 440. There seems to be some kind of threshold effect happening here that doesn't happen in Version A which I've tested out to RoO >1000 with no significant increase in solution time for Proc OptModel.

Thanks for your patience,

Gene

Version A below:

/*VERSION A*/

/* The purpose of this code is to test solution prior to implementation of solution 
provided by Rob Pratt */

/* Enter dimension of one side of Region of Optimization Matrix (RoO) centered on Yuma AA*/
%let RoO=420;
/*Define 3 possible elevation angles for each camera in radians. If a certain elevation
is to be omitted, enter a zero to indicate missing*/
options mprint symbolgen;
%LET HI=36;
%LET MID=35;
%LET LO=34;

data work.stationdataexisting;
input station $ lat long alt az elev;
datalines;
US000S 34.290983 -116.383718 .898 198.31 54.33
US000V 33.960142 -117.259044 .533 47.06 67.3
US001Q 33.960142 -117.259044 .533 129.44 68.32
US001R 34.290983 -116.383718 .898 90 45
US001E 33.656305 -117.560879 .439 135 45
;
run;
data work.stationdataexisting;
set work.stationdataexisting;
if station='US000S' then station='VS0S1R';
if station='US001R' then delete;
if station='US000V' then station='VS0V1Q';
if station='US001Q' then delete;
/*Define station grid coordinates relative to Yuma AZ*/
xpos=(long-(-114.599))*(111.321*cos(constant('pi')/180*lat));
ypos=(lat-(32.669))*111.19;
zpos= alt;
run;


data work.OrientationVectors;
	/*Define 8 possible azimuths (NO, NE, EA, SE, SO, SW, WE, NW) for each camera in radians*/
	AZN=0*constant('pi')/180;
	AZNE=45*constant('pi')/180;
	AZEA=90*constant('pi')/180;
	AZSE=135*constant('pi')/180;
	AZS=180*constant('pi')/180;
	AZSW=225*constant('pi')/180;
	AZWE=270*constant('pi')/180;
	AZNW=315*constant('pi')/180;
	ALTHI=&hi*constant('pi')/180;
	zhatHI=sin(althi);
	ALTMID=&mid*constant('pi')/180;
	zhatMID=sin(altmid);
	ALTLO=&lo*constant('pi')/180;
	zhatLO=sin(altlo);

	/*Define Unit Vectors for each camera orientation*/
	/*NORTH AZIMUTH*/
	/*AZN_ALTLO*/
	nlo_xhat=sin(azn);
	nlo_yhat=cos(azn);

	/*AZN_ALMID*/
	nmid_xhat=sin(azn);
	nmid_yhat=cos(azn);

	/*AZN_ALTHI*/
	nhi_xhat=sin(azn);
	nhi_yhat=cos(azn);

	/*NORTHEAST AZIMUTH*/
	/*AZNE_ALTLO*/
	nelo_xhat=sin(azne);
	nelo_yhat=cos(azne);

	/*AZNE_ALTMID*/
	nemid_xhat=sin(azne);
	nemid_yhat=cos(azne);

	/*AZNE_ALTHI*/
	nehi_xhat=sin(azne);
	nehi_yhat=cos(azne);
	
	/*EAST AZIMUTH*/
	/*AZEA_ALTLO*/
	ealo_xhat=sin(azea);
	ealo_yhat=cos(azea);

	/*AZEA_ALTMID*/
	eamid_xhat=sin(azea);
	eamid_yhat=cos(azea);

	/*AZEA_ALTHI*/
	eahi_xhat=sin(azea);
	eahi_yhat=cos(azea);
	
	/*SOUTHEAST AZIMUTH*/
	/*AZSE_ALTLO*/
	selo_xhat=sin(azse);
	selo_yhat=cos(azse);

	/*AZSE_ALTMID*/
	semid_xhat=sin(azse);
	semid_yhat=cos(azse);

	/*AZSE_ALTHI*/
	sehi_xhat=sin(azse);
	sehi_yhat=cos(azse);

	/*SOUTH AZIMUTH*/
	/*AZS_ALTLO*/
	slo_xhat=sin(azs);
	slo_yhat=cos(azs);

	/*AZS_ALTMID*/
	smid_xhat=sin(azs);
	smid_yhat=cos(azs);

	/*AZS_ALTHI*/
	shi_xhat=sin(azs);
	shi_yhat=cos(azs);

	/*SOUTHWEST AZIMUTH*/
	/*AZSW_ALTLO*/
	swlo_xhat=sin(azsw);
	swlo_yhat=cos(azsw);

	/*AZSW_ALTMID*/
	swmid_xhat=sin(azsw);
	swmid_yhat=cos(azsw);

	/*AZSW_ALTHI*/
	swhi_xhat=sin(azsw);
	swhi_yhat=cos(azsw);
	
	/*WEST AZIMUTH*/
	/*AZWE_ALTLO*/
	welo_xhat=sin(azwe);
	welo_yhat=cos(azwe);

	/*AZWE_ALTMID*/
	wemid_xhat=sin(azwe);
	wemid_yhat=cos(azwe);

	/*AZWE_ALTHI*/
	wehi_xhat=sin(azwe);
	wehi_yhat=cos(azwe);

	/*NORTHWEST AZIMUTH*/
	/*AZNW_ALTLO*/
	nwlo_xhat=sin(aznw);
	nwlo_yhat=cos(aznw);

	/*AZNW_ALTMID*/
	nwmid_xhat=sin(aznw);
	nwmid_yhat=cos(aznw);

	/*AZNW_ALTHI*/
	nwhi_xhat=sin(aznw);
	nwhi_yhat=cos(aznw);
run;

/* Create Coverage Matrix*/
data work.VirtualCoverageMatrix;
	format Station $8. orientation $5. gridpoint $13. k 1.;
	keep station orientation target k Tx Ty Tz;

	/*Set Target grid dimensions and grid interval*/
	Txmin=-&roo/2;
	Txmax=&roo/2;
	Tymin=-&roo/2;
	Tymax=&roo/2;
	Tzmin=70;
	Tzmax=120;
	Interval=10;

	/*Create Temporary Arrays */
	array _station(5) $20 _temporary_;
	array xyz(3,5) _temporary_;

	if _n_=1 then
		do i=1 to nobs;
			set work.stationdataexisting nobs=nobs;
			_station(i)=station;

			/*Station location*/
			xyz(1, i)=xpos;
			xyz(2, i)=ypos;
			xyz(3, i)=zpos;
		end;

	/*Temporary Array for xhat, yhat and zhat at all azimuths*/
	array vector(72) _temporary_;
	set work.orientationvectors;

	/*North Unit Vectors*/
	vector(1)=nlo_xhat;
	vector(2)=nlo_yhat;
	vector(3)=zhatlo;
	vector(4)=nmid_xhat;
	vector(5)=nmid_yhat;
	vector(6)=zhatmid;
	vector(7)=nhi_xhat;
	vector(8)=nhi_yhat;
	vector(9)=zhathi;

	/*Northeast Unit Vectors*/
	vector(10)=nelo_xhat;
	vector(11)=nelo_yhat;
	vector(12)=zhatlo;
	vector(13)=nemid_xhat;
	vector(14)=nemid_yhat;
	vector(15)=zhatmid;
	vector(16)=nehi_xhat;
	vector(17)=nehi_yhat;
	vector(18)=zhathi;
	
	/*East Unit Vectors*/
	vector(19)=ealo_xhat;
	vector(20)=ealo_yhat;
	vector(21)=zhatlo;
	vector(22)=eamid_xhat;
	vector(23)=eamid_yhat;
	vector(24)=zhatmid;
	vector(25)=eahi_xhat;
	vector(26)=eahi_yhat;
	vector(27)=zhathi;
	
	/*Southeast Unit Vectors*/
	vector(28)=selo_xhat;
	vector(29)=selo_yhat;
	vector(30)=zhatlo;
	vector(31)=semid_xhat;
	vector(32)=semid_yhat;
	vector(33)=zhatmid;
	vector(34)=sehi_xhat;
	vector(35)=sehi_yhat;
	vector(36)=zhathi;

	/*South Unit Vectors*/
	vector(37)=slo_xhat;
	vector(38)=slo_yhat;
	vector(39)=zhatlo;
	vector(40)=smid_xhat;
	vector(41)=smid_yhat;
	vector(42)=zhatmid;
	vector(43)=shi_xhat;
	vector(44)=shi_yhat;
	vector(45)=zhathi;

	/*Southwest Unit Vectors*/
	vector(46)=swlo_xhat;
	vector(47)=swlo_yhat;
	vector(48)=zhatlo;
	vector(49)=swmid_xhat;
	vector(50)=swmid_yhat;
	vector(51)=zhatmid;
	vector(52)=swhi_xhat;
	vector(53)=swhi_yhat;
	vector(54)=zhathi;
	
	/*West Unit Vectors*/
	vector(55)=welo_xhat;
	vector(56)=welo_yhat;
	vector(57)=zhatlo;
	vector(58)=wemid_xhat;
	vector(59)=wemid_yhat;
	vector(60)=zhatmid;
	vector(61)=wehi_xhat;
	vector(62)=wehi_yhat;
	vector(63)=zhathi;

	/*Northwest Unit Vectors*/
	vector(64)=nwlo_xhat;
	vector(65)=nwlo_yhat;
	vector(66)=zhatlo;
	vector(67)=nwmid_xhat;
	vector(68)=nwmid_yhat;
	vector(69)=zhatmid;
	vector(70)=nwhi_xhat;
	vector(71)=nwhi_yhat;
	vector(72)=zhathi;

	/*Start loops through Targets*/
	Target=0;

	do Tz=Tzmin to Tzmax by Interval;

		do Ty=Tymin to Tymax by Interval;

			do Tx=Txmin to Txmax by Interval;
				Target=Target + 1;
				Gridpoint=(cats(Tx, '/', Ty, '/', Tz));
				
				/* Start Loops through stations*/
				do c=1 to nobs;
					station=_station(c);

					/*station Location*/
					Cx=xyz(1, c);
					Cy=xyz(2, c);
					Cz=xyz(3, c);

					/*Create Target Vector from station to target on 3D Cartesian grid*/
					CTx=Tx-Cx;
					CTy=Ty-Cy;
					CTz=Tz-Cz;

					/*Calculate distance from station to Target in 3 Dimensions*/
					CTdist=sqrt(CTx**2+CTy**2+Ctz**2);

					/*Loop through each of 25  possible orientations of station c
					to determine whether Target(Tz, Ty, Tx) is covered (k=1) or not (k=0)*/
					do j=1 to 70 by 3;

						/* Define Position Label*/
						if j=1 then
							orientation='NOLO';

						if j=4 then
							orientation='NOMD';

						if j=7 then
							orientation='NOHI';

						if j=10 then
							orientation='NELO';

						if j=13 then
							orientation='NEMD';

						if j=16 then
							orientation='NEHI';

						if j=19 then
							orientation='EALO';

						if j=22 then
							orientation='EAMD';

						if j=25 then
							orientation='EAHI';

						if j=28 then
							orientation='SELO';

						if j=31 then
							orientation='SEMD';

						if j=34 then
							orientation='SEHI';

						if j=37 then
							orientation='SOLO';

						if j=40 then
							orientation='SOMD';

						if j=43 then
							orientation='SOHI';

						if j=46 then
							orientation='SWLO';

						if j=49 then
							orientation='SWMD';

						if j=52 then
							orientation='SWHI';
							
						if j=55 then
							orientation='WELO';

						if j=58 then
							orientation='WEMD';

						if j=61 then
							orientation='WEHI';
							
						if j=64 then
							orientation='NWLO';

						if j=67 then
							orientation='NWMD';

						if j=70 then
							orientation='NWHI';

						/*3D Unit Vector of stations at selected orientation*/
						xhat=vector(j);
						yhat=vector(j+1);
						zhat=vector(j+2);

						/*Calculate Elevation angle and  vertical FOV limits between station Location and Target Location*/
						Elev2Target=arsin(CTz/CTdist)*(180/constant('pi'));
						
						/*Fix US000S at 40 degree (zhat=sin(40)=.643)elevation and let model optimize elevation for other cameras*/
						If (find(station,'VS0S1R')) then do;
						UpperVertFoVlimit=arsin(.643)*(180/constant('pi'))+22;
						LowerVertFoVlimit=arsin(.643)*(180/constant('pi'))-22;
						end;
						else do;
						UpperVertFoVlimit=arsin(zhat)*(180/constant('pi'))+22;
						LowerVertFoVlimit=arsin(zhat)*(180/constant('pi'))-22;
						end;

						/*Test if Elevation to Target is within FoV limits*/
						If (UpperVertFoVlimit>Elev2Target and LowerVertFovlimit<Elev2Target)then
							E2T=1;
						else
							E2T=0;

						/*Calculate Angle Between station Unit Vector and Target Vector in x-y plane*/
						SepAngle2D=arcos(((xhat*CTx)+(yhat*CTy))/(sqrt(CTx**2 + CTy**2)*sqrt(xhat**2 +yhat**2)))*(180/constant('pi'));

						/* Introducing a test here to compute k for a "virtual triple camera station (VSxxx)" with a wider FoV*/
						/*0*/
						if (find(station, 'VS')) then

							/*1*/
							if ((find(station, 'VS12E')) or (find(station, 'VS89R'))) then

								/*Branch between VS stations Regular Stations*/
								/*2*/
								if (find(station, 'VS12E')) then

									/*Branch between VS12E and other VS stations*/
									/*Branch for between VS12E @ orientation NOXX and other orientations*/
									/*3*/
									if (find(orientation, 'NO')) then

										/* Compute k for VS12E @ NOXX*/
										/*4*/
										if zhat>0 and Sepangle2D LE (3*88/2) and E2T=1 and CTDist LE 320 then

											/*Sets k=1 when all conditions met*/
											do;
												k=1;
												output;
												continue;
											end;
										else
											do;

												/*Sets k=0 when all conditions not met*/
												k=0;
												output;
												continue;

												/*end 4*/
											end;

										/*else 3*/
										else
											do;

												/*Sets k=0 for other orientations of VS12E*/
												k=0;
												output;
												continue;

												/*end 3*/
											end;

										/*else 2*/
										else
											do;

												/*Computes k for VS89R by default becasue there is no other 'VS' station*/
												/*3*/
												if (find(orientation, 'SO')) then

													/*Brlanch between VS89R @ SOXX and other orientations*/
													/* Compute k for VS89R @ SOXX*/
													/*4*/
													if zhat>0 and Sepangle2D LE (3*88/2) and E2T=1 and CTDist LE 320 
														then

															/*Sets k=1 when all conditions met*/
															do;
															k=1;
															output;
															continue;
														end;
													else
														do;

															/*Sets k=0 when all conditions not met*/
															k=0;
															output;
															continue;

															/*end 4*/
														end;

													/*else 3*/
													else
														do;

															/*Sets k=0 for other orientations of VS89R*/
															k=0;
															output;
															continue;

															/*end 3*/
														end;

													/*end 2*/
												end;

											/*else 1*/
											else
												do;

													/*2*/
/* 													if ((find(orientation, 'SE')) or (find(orientation, 'NW'))) then */

														/*Branch for Other VS Stations*/
														/*3*/
														if zhat>0 and CTdist<=320 and Sepangle2D LE 2*88/2 and E2T=1 then
															do;

																/*Test for k on all orientations of other stations*/
																k=1;
																output;
																continue;
															end;

														/*else 3*/
														else
															do;
																k=0;
																output;
																continue;

																/*end 3*/
															end;

														/*else 2*/
/* 														else */
/* 															do; */
/* 																k=0; */
/* 																output; */
/* 																continue; */
/*  */
/* 																end 2 */
/* 															end; */

														/*end 1*/
													end;

												/*else 0*/
												else
													do;

														/*1*/
														if zhat>0 and CTdist<=320 and Sepangle2D LE 88/2 and E2T=1 then
															do;

																/*Test for k on all orientations of other stations*/
																k=1;
																output;
																continue;
															end;

														/*else 1*/
														else
															do;
																k=0;
																output;
																continue;

																/*end 1*/
															end;

														/*end 0*/
													end;

												/* End of loop for orientation j*/
											end;

											/*End of loop for station c*/
										end;

										/* End of loop for Tx of Targets*/
									end;

									/* End of loop for Ty of Targets*/
								end;

								/* End of loop Tz of Targets*/
							end;
							stop;
						run;



/*Subset Region of Optimization*/
data work.coveragematrix&Roo;
set work.virtualcoveragematrix;
if (abs(tx)<=&RoO/2 and abs(ty)<=&RoO/2);
run;
proc sort data=coveragematrix&RoO;
by target;
run;
/*Create a temporary datset to get a count of number of targets and to renumber target from
1 to targetcount for use by Proc Optmodel*/
data coveragematrixopt (drop=target_count);
set coveragematrix&Roo;
by target;
retain target_count 0;
if first.target then target_count=target_count+1;
target=target_count;
run;
/* Have to sort by descending target in order to get target count from obs=1*/
Proc sort data=coveragematrixopt;
By descending target;
run;
data targetcount;
set coveragematrixopt (obs=1);
keep target;
/*Re-sort by ascending target--not sure if this matters to Proc Optmodel?*/
proc sort data=coveragematrixopt;
by target;
run;
/*Use count of number of targets to create Target index file for Proc OptModel*/
data work.targets(drop=i);
set targetcount;
do i=1 to target;
target=i;
output;
end;
run;

Proc optmodel;
/* Read CoverageMatrix data into index sets*/
	set <str> stations;
	read data work.stationdataexisting  into stations=[station];
	set <str> ORIENTATIONS=/'NOHI' 'NOMD' 'NOLO' 'NEHI' 'NEMD' 'NELO' 'EAHI' 'EAMD' 'EALO' 'SEHI' 'SEMD' 'SELO' 'SOHI' 'SOMD' 'SOLO' 'SWHI' 'SWMD' 'SWLO' 'WEHI' 'WEMD' 'WELO' 'NWHI' 'NWMD' 'NWLO'/;
	set <num> TARGETS;
	read data work.targets into TARGETS=[target];
	num Matrix {stations, ORIENTATIONS, TARGETS};
	read data work.CoverageMatrixOpt into [station orientation target] Matrix=k;
	
/* Set Constants*/
	
	%Let N=5; /*number of stations*/
	num N=&N;
	%Let k=3;/* optimal number of cameras covering each target*/
	num k=&k;

/*Declare Variables*/
	/* Sadik Equation 11--Sets CHI as binary variable that is 1 when using station i and POSITION j else 0.*/ 
 
	var CHI {stationS,ORIENTATIONS} binary;
	var XIT{TARGETS} integer;
	var PSI{TARGETS} integer <=k;

/* Declare Model*/

	/*Sadik Equation 7--the function to be maximized*/
	
	max OptCoverage=sum{target in TARGETS} Psi[target]
	-sum{station in stationS, orientation in ORIENTATIONS}Chi[station,orientation];

/*Subject to Following Constraints*/

	con CHI_constraint {station in stations}: 
	sum{orientation in ORIENTATIONS}CHI[station,orientation] <=1;
	
	con XIT_CON1{target in TARGETS}:
   XIT[target]=sum{station in stations, orientation in ORIENTATIONS}
               matrix[station, orientation, target] * Chi[station,orientation];
               
	con XIT_CON2{target in TARGETS}:XIT[target]/N <= PSI[target];
	
	
	con XIT_CON3{target in TARGETS}:PSI[target] <= XIT[target];
	

/* Call Solver and Save Results*/
	solve;
	create data work.OptimalSolution_VS_&RoO from 
	[station orientation]={c in stations, o in orientations: CHI[c,o]>0.5} CHI;
	quit;
	data work.optimalsolution_VS_&RoO;
	set work.optimalsolution_VS_&RoO;
	if orientation ='NOLO' then do Az=0; Alt=&lo; end;else;
	if orientation ='NOMD' then do Az=0; Alt=&mid; end;else;
	if orientation ='NOHI' then do Az=0; Alt=&hi; end;else;
	if orientation ='NELO' then do Az=45; Alt=&lo; end;else;
	if orientation ='NEMD' then do Az=45; Alt=&mid; end;else;	
	if orientation ='NEHI' then do Az=45; Alt=&hi; end;else;
	if orientation ='EALO' then do Az=90; Alt=&lo; end;else;
	if orientation ='EAMD' then do Az=90; Alt=&mid; end;else;	
	if orientation ='EAHI' then do Az=90; Alt=&hi; end;else;
	if orientation ='SELO' then do Az=135; Alt=&lo; end;else;
	if orientation ='SEMD' then do Az=135; Alt=&mid; end;else;
	if orientation ='SEHI' then do Az=135; Alt=&hi; end;else;
	if orientation ='SOLO' then do Az=180; Alt=&lo; end;else;
	if orientation ='SOMD' then do Az=180; Alt=&mid; end;else;
	if orientation ='SOHI' then do Az=180; Alt=&hi; end;else;
	if orientation ='SWLO' then do Az=225; Alt=&lo; end;else;
	if orientation ='SWMD' then do Az=225; Alt=&mid; end;else;
	if orientation ='SWHI' then do Az=225; Alt=&hi; end;else;
	if orientation ='WELO' then do Az=270; Alt=&lo; end;else;
	if orientation ='WEMD' then do Az=270; Alt=&mid; end;else;	
	if orientation ='WEHI' then do Az=270; Alt=&hi; end;else;
	if orientation ='NWLO' then do Az=315; Alt=&lo; end;else;
	if orientation ='NWMD' then do Az=315; Alt=&mid; end;else;
	if orientation ='NWHI' then do Az=315; Alt=&hi; end;	
	run;
proc sort data=work.OptimalSolution_VS_&roo; by station;
run;

%LET Corner=%eval(&RoO/2);

/*Use Match Merge Technique to merge Stationdataexisting with OptimalSolution to create StationDataOptimal*/
/* First sort both files by Station*/
Proc sort data=work.stationdataexisting;
	by Station;
run;

Proc sort data=work.optimalsolution_vs_&RoO;
	by station;
run;

proc print data=work.optimalsolution_vs_&RoO;
Title "Optimal Solution For LO=&lo km HI=&HI km RoO=&RoO x &RoO km";
run;

Data work.stationdataoptimal socal.stationdataoptimal&RoO;
	merge work.stationdataexisting work.optimalsolution_vs_&roo;
	by Station;

	if chi=. then
		Station="US0000";
	Altrad=Alt*constant('pi')/180;
	Azrad=Az*constant('pi')/180;
	xhat=sin(azrad);
	yhat=cos(azrad);
	zhat=sin(altrad);
run;

/*Use optimized camera orientations to compute optimized FoV grid coverage*/
data work.OPTFoVGrid&roo /*optmodel.optfovgrid&RoO*/;
	keep Tx Ty Tz camera k_coverage Target;

	/*Create Temporary Arrays*/
	array _station(5) $20 _temporary_;
	array xyz(6, 5) _temporary_;

	if _n_=1 then
		do i=1 to nobs;
			set work.stationdataoptimal nobs=nobs;
			_station(i)=station;
			xyz(1, i)=xpos;
			xyz(2, i)=ypos;
			xyz(3, i)=xhat;
			xyz(4, i)=yhat;
			xyz(5, i)=zhat;
			xyz(6, i)=zpos;
		end;

	/*Set Desired FoV grid dimensions and interval*/
	Txmin=-400;
	Txmax=400;
	Tymin=-400;
	Tymax=400;
	Tzmin=70;
	Tzmax=120;
	Interval=10;

	/* Number of targets--not used elsewhere in this code at this time*/
	/*N_TARGETS=(((txmax-txmin)/interval)+1)*(((tymax-tymin)/interval)+1)*(((tzmax-tzmin)/interval)+1);*/
	/*Start loops through Target Grid*/
	Target=0;

	do Tz=Tzmin to Tzmax by Interval;

		do Ty=Tymin to Tymax by Interval;

			do Tx=Txmin to Txmax by Interval;
				Target=Target+1;

				/*Convert Tx/Ty (km relative to Big Eye intersection) to Latitude/Longitude for mapping purposes*/
/* 				Longitude=(Tx/91.09)+(-106.6293); */
/* 				Latitude=(Ty/111.19)+(35.1055); */

				/*Start loops through cameras*/
				do c=1 to nobs;
					Camera=_station(c);

					/*3D Unit Vector of Camera*/
					xhat=xyz(3, c);
					yhat=xyz(4, c);
					zhat=xyz(5, c);

					/*Camera Location*/
					Cx=xyz(1, c);
					Cy=xyz(2, c);
					Cz=xyz(6, c);

					/*Create Target Vector from Camera to target on 3D Cartesian grid*/
					CTx=Tx-Cx;
					CTy=Ty-Cy;
					CTz=Tz-Cz;

					/*Calculate distance from Camera to Target in 3 Dimensions*/
					CTdist=sqrt(CTx**2+CTy**2+Ctz**2);

					/*Test if Elevation Angle is allowed (zhat > 0)*/
					if zhat=0 then
						do;
							k_coverage=0;
							continue;
						end;
					/*Calculate Elevation angle and  vertical FOV limits between Camera Location and Target Location*/
					Elev2Target=arsin(CTz/CTdist)*(180/constant('pi'));
					
					/*Fix US000S at 40 degree (zhat=sin(40)=.643)elevation*/
						If (find(station,'VS0S1R')) then do;
						UpperVertFoVlimit=arsin(.643)*(180/constant('pi'))+22;
						LowerVertFoVlimit=arsin(.643)*(180/constant('pi'))-22;
						end;
						else do;
					UpperVertFoVlimit=arsin(xyz(5, c))*(180/constant('pi'))+22;
					LowerVertFoVlimit=arsin(xyz(5, c))*(180/constant('pi'))-22;
					end;

					/*Test if Elevation to Target is within FoV limits*/
					If (UpperVertFoVlimit>Elev2Target and LowerVertFovlimit<Elev2Target)then
						E2T=1;
					else
						E2T=0;

					/*Calculate Angle Between Camera Unit Vector and Target Vector in x-y plane*/
					SepAngle2D=arcos(((xhat*CTx)+(yhat*CTy))/(sqrt(CTx**2 + CTy**2)*sqrt(xhat**2 +yhat**2)))*(180/constant('pi'));

					/* Introducing a test here to compute k for a "virtual triple camera station (VSxxx)" with a wider FoV*/
					/*0*/
					if (find(camera, 'VS')) then

						/*1*/
						if ((find(camera, 'VS12E')) or (find(camera, 'VS89R'))) then
							if zhat>0 and Sepangle2D LE (3*88/2) and E2T=1 and CTDist LE 320 then

								/*Sets k_coverage=1 when all conditions met*/
								do;
									k_coverage=1;
									output;
									continue;
								end;
							else
								do;

									/*Sets k_coverage=0 when all conditions not met*/
									k_coverage=0;
									output;
									continue;

									/*end 4*/
								end;
							else
								do;

									/*2*/
									/*if ((find(orientation, 'SE')) or (find(orientation, 'NW'))) then*/
									/*Branch for Other VS Stations*/
									/*3*/
									if zhat>0 and CTdist<=320 and Sepangle2D LE 2*88/2 and E2T=1 then
										do;

											/*Test for k_coverage on all orientations of other stations*/
											k_coverage=1;
											output;
											continue;
										end;

									/*else 3*/
									else
										do;
											k_coverage=0;
											output;
											continue;

											/*end 3*/
										end;

									/*end 1*/
								end;

							/*else 0*/
							else
								do;

									/*1*/
									if zhat>0 and CTdist<=320 and Sepangle2D LE 88/2 and E2T=1 then
										do;

											/*Test for k_coverage on all orientations of other stations*/
											k_coverage=1;
											output;
											continue;
										end;

									/*else 1*/
									else
										do;
											k_coverage=0;
											output;
											continue;

											/*end 1*/
										end;

									/*end 0*/
								end;

							/*end; superfluous?*/
							/*End of loop for each camera for Target(n)*/
						end;

						/*End of loop for all Tx Targets in Ty-th by Tz-th layer*/
					end;

					/*End of loop for Ty-th layer*/
				end;

				/*End of loop for Tz-th level*/
			end;
			stop;
		run;

		/* Take results of FoV model and sort k_coverage by TX, TY, TZ for HeatMap Plots*/
		proc sort data=work.optfovgrid&roo;
			by tx ty tz;
		run;

		/* Take results of FoV Model and output file with sum of k_coverage by Target*/
		proc summary data=work.optfovgrid&roo;
			by tx ty tz;
			var k_coverage;
			output out=work.target_psi sum=k;
		run;

		/* Create annotation datasets per Robert Allison */
		/*Annotation Dataset for area of optimization polygon*/
		data anno_poly;
			length xc1 yc1 $10 function $50;
			linethickness=2;
			layer='front';
			x1space='datavalue';
			y1space='datavalue';
			function='polygon';
			xc1='-&Corner';
			yc1='-&Corner';
			output;
			function='polycont';
			xc1='&Corner';
			yc1='-&Corner';
			output;
			xc1='&Corner';
			yc1='&Corner';
			output;
			xc1='-&Corner';
			yc1='&Corner';
			output;

			/*Annotation Dataset for area of optimization label*/
		data anno_text;
			function='text';
			layer='front';
			anchor='left';
			label='Optimized at &RoO x &Roo km';
			x1space='datavalue';
			y1space='datavalue';
			x1=-380;
			y1=380;
			textsize=10;
			widthunit='data';
			width=400;
		run;

		/*Combine datasets for city, polygon and polygon label annotations*/
		data anno_all;
			set anno_poly anno_text socal.city_anno_yuma socal.station_yuma_anno;
		run;

		/*Build of HEATMAPPARM Code per Wicklin Blog*/
		proc format;
			value k_Fmt (Fuzz=0) /* bin k into 5 groups*/
	0="0" 1=" 1" 2="2" 3-6="3-6" 7 - high=">6";
		run;
		/*Colorizes maps according to bins defined above*/
		data order;
			length Value $11 FillColor $15;
			input raw FillColor;
			Value=put (raw, k_fmt.);
			retain ID 'SortOrder' Show 'AttrMap';
			datalines;
0 cxf8f8f8
1 verylightgray
2 verylightblue
6 moderateblue
7 moderatered
;
		run;

		/*Figures of Merit (FoM) Calculations*/
		data PSIRoO;
			set work.target_psi;

			/*Use subsetting IF to select Region of Optimization (RoO)*/
			if -&RoO/2<=tx<=&RoO/2 and -&RoO/2<=ty<=&RoO/2 and 80<=tz<=120;

			/* Limit psi so that it is always <= k per Sadik et al*/
			if k>=3 then
				psi=3;
			else
				psi=k;

			/*Compute Square of each PSI value*/
			PSI2=PSI**2;
		run;

		/*Define Formats for Proc Freq*/
		Proc format;
			value k3_6fmta Low-1='0 or 1 Camera' 2='2 Cameras' 3-6='3-6 Cameras' 
				7-high='7 or More Cameras';
		run;

		Proc freq data=psiroo;
			title height=12pt "Distribution of Camera Coverage";
			tables k/;
			format k k3_6fmta.;
			label k="Over the Region of Optimization:" &Roo "KM by " &RoO "KM" @ &lo &mid &hi;
		run;

		Proc format;
			value k3_6fmtb 1='1 Camera' 2='2 Cameras' 3-6='3-6 Cameras' 
				7-high='7 or More Cameras';

		Proc freq data=target_psi;
			where k GE 1;
			title height=12pt "Distribution of Camera Coverage";
			tables k/;
			format k k3_6fmtb.;
			label k="Over the Region of Coverage";
		run;

		/*Fairness and Balanced Index from Sadik et al.*/
		/* Compute sum of PSI and PSI**2*/
		proc summary data=psiroo;
			var psi psi2;
			output out=work.psisums sum=psisum psi2sum;
		run;

		/*Compute Fairness Index and Balanced Index*/
		data _null_;
			file print;
			set psisums;

			/*Jains Fairness Index for all targets*/
			num=psisum**2;
			den=_freq_*psi2sum;
			FI=num/den;

			/*Coverage Index: Fraction of total attainable coverage for k=3*/
			CI=psisum/(_freq_*3);

			/*Sadik at al Balancing Index Calculation*/
			BI=FI*CI;

			/*  Maximize Product of  Volume of Sky in RoO and BI = Area Weighted Balancing Index?*/
			AWBI=_freq_*BI;
			put _freq_=;
			title height=12 pt 
				'Fairness, Coverage and Balancing Indexes per Sadik et al.';
			put FI=bestd9.;
			put CI=bestd9.;
			put BI=bestd9.;
			put AWBI=bestd9.;
		run;

		ods graphics / reset width=6.4in height=6.4in imagemap;
				
		
		proc sgplot data=work.target_psi dattrmap=order sganno=anno_all aspect=1;
			where tz=100;
			title height=12pt "Grouped Camera Coverage at 100km";
		
			title2 height=10pt "LO=&lo MID=&mid HI=&HI";
			format k k_fmt.;
			heatmapparm x=tx y=ty colorgroup=k/ attrid=SortOrder;
			keylegend / title="Cameras Viewing Each Cell";
			xaxis label='KM' valuesrotate=vertical grid values=(-400 to 400 by 25) 
				fitpolicy=rotatethin;
			yaxis label='KM' grid values=(-400 to 400 by 25) fitpolicy=thin;
		run;
	V

Version B below:

/*VERSION B*/

/* The purpose of this code is to test solution provided by Rob Pratt to optimize 3-camera coverage*/

/* Enter dimension of one side of Region of Optimization Matrix (RoO) centered on Yuma AZ*/
options mprint symbolgen;
%let RoO=420;

/*Define 3 possible elevation angles for each camera in radians. If a certain elevation
is to be omitted, enter a zero to indicate missing*/
%LET HI=36;
%LET MID=35;
%LET LO=34;

data work.stationdataexisting;
input station $ lat long alt az elev;
datalines;
US000S 34.290983 -116.383718 .898 198.31 54.33
US000V 33.960142 -117.259044 .533 47.06 67.3
US001Q 33.960142 -117.259044 .533 129.44 68.32
US001R 34.290983 -116.383718 .898 90 45
US001E 33.656305 -117.560879 .439 135 45
;
run;
data work.stationdataexisting;
set work.stationdataexisting;
if station='US000S' then station='VS0S1R';
if station='US001R' then delete;
if station='US000V' then station='VS0V1Q';
if station='US001Q' then delete;
/*Define station grid coordinates relative to Yuma AZ*/
xpos=(long-(-114.599))*(111.321*cos(constant('pi')/180*lat));
ypos=(lat-(32.669))*111.19;
zpos= alt;
run;


data work.OrientationVectors;
	/*Define 8 possible azimuths (NO, NE, EA, SE, SO, SW, WE, NW) for each camera in radians*/
	AZN=0*constant('pi')/180;
	AZNE=45*constant('pi')/180;
	AZEA=90*constant('pi')/180;
	AZSE=135*constant('pi')/180;
	AZS=180*constant('pi')/180;
	AZSW=225*constant('pi')/180;
	AZWE=270*constant('pi')/180;
	AZNW=315*constant('pi')/180;
	ALTHI=&hi*constant('pi')/180;
	zhatHI=sin(althi);
	ALTMID=&mid*constant('pi')/180;
	zhatMID=sin(altmid);
	ALTLO=&lo*constant('pi')/180;
	zhatLO=sin(altlo);

	/*Define Unit Vectors for each camera orientation*/
	/*NORTH AZIMUTH*/
	/*AZN_ALTLO*/
	nlo_xhat=sin(azn);
	nlo_yhat=cos(azn);

	/*AZN_ALMID*/
	nmid_xhat=sin(azn);
	nmid_yhat=cos(azn);

	/*AZN_ALTHI*/
	nhi_xhat=sin(azn);
	nhi_yhat=cos(azn);

	/*NORTHEAST AZIMUTH*/
	/*AZNE_ALTLO*/
	nelo_xhat=sin(azne);
	nelo_yhat=cos(azne);

	/*AZNE_ALTMID*/
	nemid_xhat=sin(azne);
	nemid_yhat=cos(azne);

	/*AZNE_ALTHI*/
	nehi_xhat=sin(azne);
	nehi_yhat=cos(azne);
	
	/*EAST AZIMUTH*/
	/*AZEA_ALTLO*/
	ealo_xhat=sin(azea);
	ealo_yhat=cos(azea);

	/*AZEA_ALTMID*/
	eamid_xhat=sin(azea);
	eamid_yhat=cos(azea);

	/*AZEA_ALTHI*/
	eahi_xhat=sin(azea);
	eahi_yhat=cos(azea);
	
	/*SOUTHEAST AZIMUTH*/
	/*AZSE_ALTLO*/
	selo_xhat=sin(azse);
	selo_yhat=cos(azse);

	/*AZSE_ALTMID*/
	semid_xhat=sin(azse);
	semid_yhat=cos(azse);

	/*AZSE_ALTHI*/
	sehi_xhat=sin(azse);
	sehi_yhat=cos(azse);

	/*SOUTH AZIMUTH*/
	/*AZS_ALTLO*/
	slo_xhat=sin(azs);
	slo_yhat=cos(azs);

	/*AZS_ALTMID*/
	smid_xhat=sin(azs);
	smid_yhat=cos(azs);

	/*AZS_ALTHI*/
	shi_xhat=sin(azs);
	shi_yhat=cos(azs);

	/*SOUTHWEST AZIMUTH*/
	/*AZSW_ALTLO*/
	swlo_xhat=sin(azsw);
	swlo_yhat=cos(azsw);

	/*AZSW_ALTMID*/
	swmid_xhat=sin(azsw);
	swmid_yhat=cos(azsw);

	/*AZSW_ALTHI*/
	swhi_xhat=sin(azsw);
	swhi_yhat=cos(azsw);
	
	/*WEST AZIMUTH*/
	/*AZWE_ALTLO*/
	welo_xhat=sin(azwe);
	welo_yhat=cos(azwe);

	/*AZWE_ALTMID*/
	wemid_xhat=sin(azwe);
	wemid_yhat=cos(azwe);

	/*AZWE_ALTHI*/
	wehi_xhat=sin(azwe);
	wehi_yhat=cos(azwe);

	/*NORTHWEST AZIMUTH*/
	/*AZNW_ALTLO*/
	nwlo_xhat=sin(aznw);
	nwlo_yhat=cos(aznw);

	/*AZNW_ALTMID*/
	nwmid_xhat=sin(aznw);
	nwmid_yhat=cos(aznw);

	/*AZNW_ALTHI*/
	nwhi_xhat=sin(aznw);
	nwhi_yhat=cos(aznw);
run;

/* Create Coverage Matrix*/
data work.VirtualCoverageMatrix;
	format Station $8. orientation $5. gridpoint $13. k 1.;
	keep station orientation target k Tx Ty Tz;

	/*Set Target grid dimensions and grid interval*/
	Txmin=-&roo/2;
	Txmax=&roo/2;
	Tymin=-&roo/2;
	Tymax=&roo/2;
	Tzmin=70;
	Tzmax=120;
	Interval=10;

	/*Create Temporary Arrays */
	array _station(5) $20 _temporary_;
	array xyz(3,5) _temporary_;

	if _n_=1 then
		do i=1 to nobs;
			set work.stationdataexisting nobs=nobs;
			_station(i)=station;

			/*Station location*/
			xyz(1, i)=xpos;
			xyz(2, i)=ypos;
			xyz(3, i)=zpos;
		end;

	/*Temporary Array for xhat, yhat and zhat at all azimuths*/
	array vector(72) _temporary_;
	set work.orientationvectors;

	/*North Unit Vectors*/
	vector(1)=nlo_xhat;
	vector(2)=nlo_yhat;
	vector(3)=zhatlo;
	vector(4)=nmid_xhat;
	vector(5)=nmid_yhat;
	vector(6)=zhatmid;
	vector(7)=nhi_xhat;
	vector(8)=nhi_yhat;
	vector(9)=zhathi;

	/*Northeast Unit Vectors*/
	vector(10)=nelo_xhat;
	vector(11)=nelo_yhat;
	vector(12)=zhatlo;
	vector(13)=nemid_xhat;
	vector(14)=nemid_yhat;
	vector(15)=zhatmid;
	vector(16)=nehi_xhat;
	vector(17)=nehi_yhat;
	vector(18)=zhathi;
	
	/*East Unit Vectors*/
	vector(19)=ealo_xhat;
	vector(20)=ealo_yhat;
	vector(21)=zhatlo;
	vector(22)=eamid_xhat;
	vector(23)=eamid_yhat;
	vector(24)=zhatmid;
	vector(25)=eahi_xhat;
	vector(26)=eahi_yhat;
	vector(27)=zhathi;
	
	/*Southeast Unit Vectors*/
	vector(28)=selo_xhat;
	vector(29)=selo_yhat;
	vector(30)=zhatlo;
	vector(31)=semid_xhat;
	vector(32)=semid_yhat;
	vector(33)=zhatmid;
	vector(34)=sehi_xhat;
	vector(35)=sehi_yhat;
	vector(36)=zhathi;

	/*South Unit Vectors*/
	vector(37)=slo_xhat;
	vector(38)=slo_yhat;
	vector(39)=zhatlo;
	vector(40)=smid_xhat;
	vector(41)=smid_yhat;
	vector(42)=zhatmid;
	vector(43)=shi_xhat;
	vector(44)=shi_yhat;
	vector(45)=zhathi;

	/*Southwest Unit Vectors*/
	vector(46)=swlo_xhat;
	vector(47)=swlo_yhat;
	vector(48)=zhatlo;
	vector(49)=swmid_xhat;
	vector(50)=swmid_yhat;
	vector(51)=zhatmid;
	vector(52)=swhi_xhat;
	vector(53)=swhi_yhat;
	vector(54)=zhathi;
	
	/*West Unit Vectors*/
	vector(55)=welo_xhat;
	vector(56)=welo_yhat;
	vector(57)=zhatlo;
	vector(58)=wemid_xhat;
	vector(59)=wemid_yhat;
	vector(60)=zhatmid;
	vector(61)=wehi_xhat;
	vector(62)=wehi_yhat;
	vector(63)=zhathi;

	/*Northwest Unit Vectors*/
	vector(64)=nwlo_xhat;
	vector(65)=nwlo_yhat;
	vector(66)=zhatlo;
	vector(67)=nwmid_xhat;
	vector(68)=nwmid_yhat;
	vector(69)=zhatmid;
	vector(70)=nwhi_xhat;
	vector(71)=nwhi_yhat;
	vector(72)=zhathi;

	/*Start loops through Targets*/
	Target=0;

	do Tz=Tzmin to Tzmax by Interval;

		do Ty=Tymin to Tymax by Interval;

			do Tx=Txmin to Txmax by Interval;
				Target=Target + 1;
				Gridpoint=(cats(Tx, '/', Ty, '/', Tz));

				/* Start Loops through stations*/
				do c=1 to nobs;
					station=_station(c);

					/*station Location*/
					Cx=xyz(1, c);
					Cy=xyz(2, c);
					Cz=xyz(3, c);

					/*Create Target Vector from station to target on 3D Cartesian grid*/
					CTx=Tx-Cx;
					CTy=Ty-Cy;
					CTz=Tz-Cz;

					/*Calculate distance from station to Target in 3 Dimensions*/
					CTdist=sqrt(CTx**2+CTy**2+Ctz**2);

					/*Loop through each of 24  possible orientations of station c
					to determine whether Target(Tz, Ty, Tx) is covered (k=1) or not (k=0)*/
					do j=1 to 70 by 3;

						/* Define Position Label*/
						if j=1 then
							orientation='NOLO';

						if j=4 then
							orientation='NOMD';

						if j=7 then
							orientation='NOHI';

						if j=10 then
							orientation='NELO';

						if j=13 then
							orientation='NEMD';

						if j=16 then
							orientation='NEHI';

						if j=19 then
							orientation='EALO';

						if j=22 then
							orientation='EAMD';

						if j=25 then
							orientation='EAHI';

						if j=28 then
							orientation='SELO';

						if j=31 then
							orientation='SEMD';

						if j=34 then
							orientation='SEHI';

						if j=37 then
							orientation='SOLO';

						if j=40 then
							orientation='SOMD';

						if j=43 then
							orientation='SOHI';

						if j=46 then
							orientation='SWLO';

						if j=49 then
							orientation='SWMD';

						if j=52 then
							orientation='SWHI';
							
						if j=55 then
							orientation='WELO';

						if j=58 then
							orientation='WEMD';

						if j=61 then
							orientation='WEHI';
							
						if j=64 then
							orientation='NWLO';

						if j=67 then
							orientation='NWMD';

						if j=70 then
							orientation='NWHI';

						/*3D Unit Vector of stations at selected orientation*/
						xhat=vector(j);
						yhat=vector(j+1);
						zhat=vector(j+2);

						/*Test if Elevation Angle is allowed (zhat > 0)*/
						/*if zhat = 0 then do;
						k=0;
						continue;
						end;*/
						/*Calculate Elevation angle and  vertical FOV limits between station Location and Target Location*/
						Elev2Target=arsin(CTz/CTdist)*(180/constant('pi'));
						
						/*Fix US000S at 40 degree (zhat=sin(40)=.643)elevation and let model optimize elevation for other cameras*/
						If (find(station,'VS0S1R')) then do;
						UpperVertFoVlimit=arsin(.643)*(180/constant('pi'))+22;
						LowerVertFoVlimit=arsin(.643)*(180/constant('pi'))-22;
						end;
						else do;
						UpperVertFoVlimit=arsin(zhat)*(180/constant('pi'))+22;
						LowerVertFoVlimit=arsin(zhat)*(180/constant('pi'))-22;
						end;

						/*Test if Elevation to Target is within FoV limits*/
						If (UpperVertFoVlimit>Elev2Target and LowerVertFovlimit<Elev2Target)then
							E2T=1;
						else
							E2T=0;

						/*Calculate Angle Between station Unit Vector and Target Vector in x-y plane*/
						SepAngle2D=arcos(((xhat*CTx)+(yhat*CTy))/(sqrt(CTx**2 + CTy**2)*sqrt(xhat**2 +yhat**2)))*(180/constant('pi'));

						/* Introducing a test here to compute k for a "virtual triple camera station (VSxxx)" with a wider FoV*/
						/*0*/
						if (find(station, 'VS')) then

							/*1*/
							if ((find(station, 'VS12E')) or (find(station, 'VS89R'))) then

								/*Branch between VS stations Regular Stations*/
								/*2*/
								if (find(station, 'VS12E')) then

									/*Branch between VS12E and other VS stations*/
									/*Branch for between VS12E @ orientation NOXX and other orientations*/
									/*3*/
									if (find(orientation, 'NO')) then

										/* Compute k for VS12E @ NOXX*/
										/*4*/
										if zhat>0 and Sepangle2D LE (3*88/2) and E2T=1 and CTDist LE 320 then

											/*Sets k=1 when all conditions met*/
											do;
												k=1;
												output;
												continue;
											end;
										else
											do;

												/*Sets k=0 when all conditions not met*/
												k=0;
												output;
												continue;

												/*end 4*/
											end;

										/*else 3*/
										else
											do;

												/*Sets k=0 for other orientations of VS12E*/
												k=0;
												output;
												continue;

												/*end 3*/
											end;

										/*else 2*/
										else
											do;

												/*Computes k for VS89R by default becasue there is no other 'VS' station*/
												/*3*/
												if (find(orientation, 'SO')) then

													/*Brlanch between VS89R @ SOXX and other orientations*/
													/* Compute k for VS89R @ SOXX*/
													/*4*/
													if zhat>0 and Sepangle2D LE (3*88/2) and E2T=1 and CTDist LE 320 
														then

															/*Sets k=1 when all conditions met*/
															do;
															k=1;
															output;
															continue;
														end;
													else
														do;

															/*Sets k=0 when all conditions not met*/
															k=0;
															output;
															continue;

															/*end 4*/
														end;

													/*else 3*/
													else
														do;

															/*Sets k=0 for other orientations of VS89R*/
															k=0;
															output;
															continue;

															/*end 3*/
														end;

													/*end 2*/
												end;

											/*else 1*/
											else
												do;

													/*2*/
/* 													if ((find(orientation, 'SE')) or (find(orientation, 'NW'))) then */

														/*Branch for Other VS Stations*/
														/*3*/
														if zhat>0 and CTdist<=320 and Sepangle2D LE 2*88/2 and E2T=1 then
															do;

																/*Test for k on all orientations of other stations*/
																k=1;
																output;
																continue;
															end;

														/*else 3*/
														else
															do;
																k=0;
																output;
																continue;

																/*end 3*/
															end;

														/*else 2*/
/* 														else */
/* 															do; */
/* 																k=0; */
/* 																output; */
/* 																continue; */
/*  */
/* 																end 2 */
/* 															end; */

														/*end 1*/
													end;

												/*else 0*/
												else
													do;

														/*1*/
														if zhat>0 and CTdist<=320 and Sepangle2D LE 88/2 and E2T=1 then
															do;

																/*Test for k on all orientations of other stations*/
																k=1;
																output;
																continue;
															end;

														/*else 1*/
														else
															do;
																k=0;
																output;
																continue;

																/*end 1*/
															end;

														/*end 0*/
													end;

												/* End of loop for orientation j*/
											end;

											/*End of loop for station c*/
										end;

										/* End of loop for Tx of Targets*/
									end;

									/* End of loop for Ty of Targets*/
								end;

								/* End of loop Tz of Targets*/
							end;
							stop;
						run;

/*Subset Region of Optimization*/
data work.coveragematrix&Roo;
set work.virtualcoveragematrix;
if (abs(tx)<=&RoO/2 and abs(ty)<=&RoO/2);
run;
proc sort data=coveragematrix&RoO;
by target;
run;
/*Create a temporary datset to get a count of number of targets and to renumber target from
1 to targetcount for use by Proc Optmodel*/
data coveragematrixopt (drop=target_count);
set coveragematrix&Roo;
by target;
retain target_count 0;
if first.target then target_count=target_count+1;
target=target_count;
run;
/* Have to sort by descending target in order to get target count from obs=1*/
Proc sort data=coveragematrixopt;
By descending target;
run;
data targetcount;
set coveragematrixopt (obs=1);
keep target;
/*Re-sort by ascending target--not sure if this matters to Proc Optmodel?*/
proc sort data=coveragematrixopt;
by target;
run;
/*Use count of number of targets to create Target index file for Proc OptModel*/
data work.targets(drop=i);
set targetcount;
do i=1 to target;
target=i;
output;
end;
run;

Proc optmodel;
/* Read CoverageMatrix data into index sets*/
	set <str> stations;
	read data work.stationdataexisting  into stations=[station];
	set <str> ORIENTATIONS=/'NOHI' 'NOMD' 'NOLO' 'NEHI' 'NEMD' 'NELO' 'EAHI' 'EAMD' 'EALO' 'SEHI' 'SEMD' 'SELO' 'SOHI' 'SOMD' 'SOLO' 'SWHI' 'SWMD' 'SWLO' 'WEHI' 'WEMD' 'WELO' 'NWHI' 'NWMD' 'NWLO'/;
	set <num> TARGETS;
	read data work.targets into TARGETS=[target];
	num Matrix {stations, ORIENTATIONS, TARGETS};
	read data work.CoverageMatrixOpt into [station orientation target] Matrix=k;
	
/* Set Constants*/
	
	%Let N=5; /*number of stations*/
	num N=&N;
	%Let k=3;/* optimal number of cameras covering each target*/
	num k=&k;

/*Declare Variables*/
 
	var CHI {stationS,ORIENTATIONS} binary;
	var IsTriplyCovered {TARGETS} binary;
	var XIT{TARGETS} integer;

	/* Declare Model*/
	
	/*Rob Pratt's suggested function to be maximized*/

	max NumTriplyCovered = sum{t in TARGETS} IsTriplyCovered[t];


/*Subject to Following Constraints*/
	con CHI_constraint {station in stations}: 
	sum{orientation in ORIENTATIONS}CHI[station, orientation] <=1;
	
	con TriplyCoveredCon {t in TARGETS}: XIT[t] >=3 * IsTriplyCovered[t];
	
	con XIT_CON1{target in TARGETS}:
   XIT[target]=sum{station in stations, orientation in ORIENTATIONS} 
		matrix[station, orientation, target] * Chi[station, orientation];
	
	/* Call Solver and Save Results*/

	solve;
	create data work.OptimalSolution_VS_&RoO from 
	[station orientation]={c in stations, o in orientations: CHI[c,o]>0.5} CHI;
	quit;
	data work.optimalsolution_VS_&RoO;
	set work.optimalsolution_VS_&RoO;
	if orientation ='NOLO' then do Az=0; Alt=&lo; end;else;
	if orientation ='NOMD' then do Az=0; Alt=&mid; end;else;
	if orientation ='NOHI' then do Az=0; Alt=&hi; end;else;
	if orientation ='NELO' then do Az=45; Alt=&lo; end;else;
	if orientation ='NEMD' then do Az=45; Alt=&mid; end;else;	
	if orientation ='NEHI' then do Az=45; Alt=&hi; end;else;
	if orientation ='EALO' then do Az=90; Alt=&lo; end;else;
	if orientation ='EAMD' then do Az=90; Alt=&mid; end;else;	
	if orientation ='EAHI' then do Az=90; Alt=&hi; end;else;
	if orientation ='SELO' then do Az=135; Alt=&lo; end;else;
	if orientation ='SEMD' then do Az=135; Alt=&mid; end;else;
	if orientation ='SEHI' then do Az=135; Alt=&hi; end;else;
	if orientation ='SOLO' then do Az=180; Alt=&lo; end;else;
	if orientation ='SOMD' then do Az=180; Alt=&mid; end;else;
	if orientation ='SOHI' then do Az=180; Alt=&hi; end;else;
	if orientation ='SWLO' then do Az=225; Alt=&lo; end;else;
	if orientation ='SWMD' then do Az=225; Alt=&mid; end;else;
	if orientation ='SWHI' then do Az=225; Alt=&hi; end;else;
	if orientation ='WELO' then do Az=270; Alt=&lo; end;else;
	if orientation ='WEMD' then do Az=270; Alt=&mid; end;else;	
	if orientation ='WEHI' then do Az=270; Alt=&hi; end;else;
	if orientation ='NWLO' then do Az=315; Alt=&lo; end;else;
	if orientation ='NWMD' then do Az=315; Alt=&mid; end;else;
	if orientation ='NWHI' then do Az=315; Alt=&hi; end;	
	run;
proc sort data=work.OptimalSolution_VS_&roo; by station;
run;

%LET Corner=%eval(&RoO/2);

Proc sort data=work.stationdataexisting;
	by Station;
run;

Proc sort data=work.optimalsolution_vs_&RoO;
	by station;
run;

proc print data=work.optimalsolution_vs_&RoO;
Title "Optimal Solution For LO=&lo km HI=&HI km RoO=&RoO x &RoO km";
run;

Data work.stationdataoptimal socal.stationdataoptimal&RoO;
	merge work.stationdataexisting work.optimalsolution_vs_&roo;
	by Station;

	if chi=. then
		Station="US0000";
	Altrad=Alt*constant('pi')/180;
	Azrad=Az*constant('pi')/180;
	xhat=sin(azrad);
	yhat=cos(azrad);
	zhat=sin(altrad);
run;

/*Use optimized camera orientations to compute optimized FoV grid coverage*/
data work.OPTFoVGrid&roo /*optmodel.optfovgrid&RoO*/;
	keep Tx Ty Tz camera k_coverage Target;

	/*Create Temporary Arrays*/
	array _station(5) $20 _temporary_;
	array xyz(6, 5) _temporary_;

	if _n_=1 then
		do i=1 to nobs;
			set work.stationdataoptimal nobs=nobs;
			_station(i)=station;
			xyz(1, i)=xpos;
			xyz(2, i)=ypos;
			xyz(3, i)=xhat;
			xyz(4, i)=yhat;
			xyz(5, i)=zhat;
			xyz(6, i)=zpos;
		end;

	/*Set Desired FoV grid dimensions and interval*/
	Txmin=-400;
	Txmax=400;
	Tymin=-400;
	Tymax=400;
	Tzmin=70;
	Tzmax=120;
	Interval=10;
	Target=0;

	do Tz=Tzmin to Tzmax by Interval;

		do Ty=Tymin to Tymax by Interval;

			do Tx=Txmin to Txmax by Interval;
				Target=Target+1;

				/*Start loops through cameras*/
				do c=1 to nobs;
					Camera=_station(c);

					/*3D Unit Vector of Camera*/
					xhat=xyz(3, c);
					yhat=xyz(4, c);
					zhat=xyz(5, c);

					/*Camera Location*/
					Cx=xyz(1, c);
					Cy=xyz(2, c);
					Cz=xyz(6, c);

					/*Create Target Vector from Camera to target on 3D Cartesian grid*/
					CTx=Tx-Cx;
					CTy=Ty-Cy;
					CTz=Tz-Cz;

					/*Calculate distance from Camera to Target in 3 Dimensions*/
					CTdist=sqrt(CTx**2+CTy**2+Ctz**2);

					/*Test if Elevation Angle is allowed (zhat > 0)*/
					if zhat=0 then
						do;
							k_coverage=0;
							continue;
						end;
					/*Calculate Elevation angle and  vertical FOV limits between Camera Location and Target Location*/
					Elev2Target=arsin(CTz/CTdist)*(180/constant('pi'));
					
					/*Fix US000S at 40 degree (zhat=sin(40)=.643)elevation*/
						If (find(station,'VS0S1R')) then do;
						UpperVertFoVlimit=arsin(.643)*(180/constant('pi'))+22;
						LowerVertFoVlimit=arsin(.643)*(180/constant('pi'))-22;
						end;
						else do;
					UpperVertFoVlimit=arsin(xyz(5, c))*(180/constant('pi'))+22;
					LowerVertFoVlimit=arsin(xyz(5, c))*(180/constant('pi'))-22;
					end;

					/*Test if Elevation to Target is within FoV limits*/
					If (UpperVertFoVlimit>Elev2Target and LowerVertFovlimit<Elev2Target)then
						E2T=1;
					else
						E2T=0;

					/*Calculate Angle Between Camera Unit Vector and Target Vector in x-y plane*/
					SepAngle2D=arcos(((xhat*CTx)+(yhat*CTy))/(sqrt(CTx**2 + CTy**2)*sqrt(xhat**2 +yhat**2)))*(180/constant('pi'));

					/* Introducing a test here to compute k for a "virtual triple camera station (VSxxx)" with a wider FoV*/
					/*0*/
					if (find(camera, 'VS')) then

						/*1*/
						if ((find(camera, 'VS12E')) or (find(camera, 'VS89R'))) then
							if zhat>0 and Sepangle2D LE (3*88/2) and E2T=1 and CTDist LE 320 then

								/*Sets k_coverage=1 when all conditions met*/
								do;
									k_coverage=1;
									output;
									continue;
								end;
							else
								do;

									/*Sets k_coverage=0 when all conditions not met*/
									k_coverage=0;
									output;
									continue;

									/*end 4*/
								end;
							else
								do;

									/*2*/
									/*if ((find(orientation, 'SE')) or (find(orientation, 'NW'))) then*/
									/*Branch for Other VS Stations*/
									/*3*/
									if zhat>0 and CTdist<=320 and Sepangle2D LE 2*88/2 and E2T=1 then
										do;

											/*Test for k_coverage on all orientations of other stations*/
											k_coverage=1;
											output;
											continue;
										end;

									/*else 3*/
									else
										do;
											k_coverage=0;
											output;
											continue;

											/*end 3*/
										end;

									/*end 1*/
								end;

							/*else 0*/
							else
								do;

									/*1*/
									if zhat>0 and CTdist<=320 and Sepangle2D LE 88/2 and E2T=1 then
										do;

											/*Test for k_coverage on all orientations of other stations*/
											k_coverage=1;
											output;
											continue;
										end;

									/*else 1*/
									else
										do;
											k_coverage=0;
											output;
											continue;

											/*end 1*/
										end;

									/*end 0*/
								end;

							/*end; superfluous?*/
							/*End of loop for each camera for Target(n)*/
						end;

						/*End of loop for all Tx Targets in Ty-th by Tz-th layer*/
					end;

					/*End of loop for Ty-th layer*/
				end;

				/*End of loop for Tz-th level*/
			end;
			stop;
		run;

		/* Take results of FoV model and sort k_coverage by TX, TY, TZ for HeatMap Plots*/
		proc sort data=work.optfovgrid&roo;
			by tx ty tz;
		run;

		/* Take results of FoV Model and output file with sum of k_coverage by Target*/
		proc summary data=work.optfovgrid&roo;
			by tx ty tz;
			var k_coverage;
			output out=work.target_psi sum=k;
		run;

		/*Annotation Dataset for area of optimization polygon*/
		data anno_poly;
			length xc1 yc1 $10 function $50;
			linethickness=2;
			layer='front';
			x1space='datavalue';
			y1space='datavalue';
			function='polygon';
			xc1='-&Corner';
			yc1='-&Corner';
			output;
			function='polycont';
			xc1='&Corner';
			yc1='-&Corner';
			output;
			xc1='&Corner';
			yc1='&Corner';
			output;
			xc1='-&Corner';
			yc1='&Corner';
			output;

			/*Annotation Dataset for area of optimization label*/
		data anno_text;
			function='text';
			layer='front';
			anchor='left';
			label='Optimized at &RoO x &Roo km';
			x1space='datavalue';
			y1space='datavalue';
			x1=-380;
			y1=380;
			textsize=10;
			widthunit='data';
			width=400;
		run;

		/*Combine datasets for city, polygon and polygon label annotations*/
		data anno_all;
			set anno_poly anno_text socal.city_anno_yuma socal.station_yuma_anno;
		run;

		/*Build of HEATMAPPARM Code per Wicklin Blog*/
		proc format;
			value k_Fmt (Fuzz=0) /* bin k into 5 groups*/
	0="0" 1=" 1" 2="2" 3-6="3-6" 7 - high=">6";
		run;

		data order;
			length Value $11 FillColor $15;
			input raw FillColor;
			Value=put (raw, k_fmt.);
			retain ID 'SortOrder' Show 'AttrMap';
			datalines;
0 cxf8f8f8
1 verylightgray
2 verylightblue
6 moderateblue
7 moderatered
;
		run;

		/*Figures of Merit (FoM) Calculations*/
		data PSIRoO;
			set work.target_psi;

			/*Use subsetting IF to select Region of Optimization (RoO)*/
			if -&RoO/2<=tx<=&RoO/2 and -&RoO/2<=ty<=&RoO/2 and 80<=tz<=120;

			/* Limit psi so that it is always <= k per Sadik et al*/
			if k>=3 then
				psi=3;
			else
				psi=k;

			/*Compute Square of each PSI value*/
			PSI2=PSI**2;
		run;

		/*Define Formats for Proc Freq*/
		Proc format;
			value k3_6fmta Low-1='0 or 1 Camera' 2='2 Cameras' 3-6='3-6 Cameras' 
				7-high='7 or More Cameras';
		run;

		Proc freq data=psiroo;
			title height=12pt "Distribution of Camera Coverage";
			tables k/;
			format k k3_6fmta.;
			label k="Over the Region of Optimization:" &Roo "KM by " &RoO "KM" @ &lo &mid &hi;
		run;

		Proc format;
			value k3_6fmtb 1='1 Camera' 2='2 Cameras' 3-6='3-6 Cameras' 
				7-high='7 or More Cameras';

		Proc freq data=target_psi;
			where k GE 1;
			title height=12pt "Distribution of Camera Coverage";
			tables k/;
			format k k3_6fmtb.;
			label k="Over the Region of Coverage";
		run;

		/*Fairness and Balanced Index from Sadik et al.*/
		/* Compute sum of PSI and PSI**2*/
		proc summary data=psiroo;
			var psi psi2;
			output out=work.psisums sum=psisum psi2sum;
		run;

		/*Compute Fairness Index and Balanced Index*/
		data _null_;
			file print;
			set psisums;

			/*Jains Fairness Index for all targets*/
			num=psisum**2;
			den=_freq_*psi2sum;
			FI=num/den;

			/*Coverage Index: Fraction of total attainable coverage for k=3*/
			CI=psisum/(_freq_*3);

			/*Sadik at al Balancing Index Calculation*/
			BI=FI*CI;

			/*  Maximize Product of  Volume of Sky in RoO and BI = Area Weighted Balancing Index?*/
			AWBI=_freq_*BI;
			put _freq_=;
			title height=12 pt 
				'Fairness, Coverage and Balancing Indexes per Sadik et al.';
			put FI=bestd9.;
			put CI=bestd9.;
			put BI=bestd9.;
			put AWBI=bestd9.;
		run;

		ods graphics / reset width=6.4in height=6.4in imagemap;
		
		proc sgplot data=work.target_psi dattrmap=order sganno=anno_all aspect=1;
			where tz=100;
			title height=12pt "Grouped Camera Coverage at 100km";
		
			title2 height=10pt "LO=&lo MID=&mid HI=&HI";
			format k k_fmt.;
			heatmapparm x=tx y=ty colorgroup=k/ attrid=SortOrder;
			keylegend / title="Cameras Viewing Each Cell";
			xaxis label='KM' valuesrotate=vertical grid values=(-400 to 400 by 25) 
				fitpolicy=rotatethin;
			yaxis label='KM' grid values=(-400 to 400 by 25) fitpolicy=thin;
		run;

RobPratt · Posted 06-21-2021 05:13 PM

For Version B with 420, the MILP presolver is able to remove almost all of the variables and constraints. For Version B with 440, it is not able to do so, for some reason that I cannot yet explain. You can save the presolver some work in this case by replacing the explicit variable XIT and constraint XIT_CON1 with the following implicit variable:

	impvar XIT {target in TARGETS} =
      sum{station in stations, orientation in ORIENTATIONS} 
   		matrix[station, orientation, target] * Chi[station, orientation];

Remember to declare the TriplyCoveredCon constraint after this IMPVAR statement.

With this change for Version B, both 420 and 440 (and even larger values of RoO) solve instantly.

genemroz · Posted 06-21-2021 11:11 PM

Rob,

Thanks for taking the time to work this vexing issue. If you learn why the MILB presolver didn't work as expected, please let me know as now I'm curious. Your solution to use an implicit variable appears to work well and I will mark it as 'accepted'. I will be putting it through additional tests (increasing the number of cameras) in the coming days. If I run into any issues, I'll post again on this community board.

Thanks again,

Gene

pchristophel · Posted 06-23-2021 12:24 PM

Hello Gene,

I looked into the inner workings to make sure there is nothing wrong. As it turns out the presolver tries to not spend excessive time substituting out fairly long columns (more than 5000 non-zeros) if it can't be sure that this will be faster than just starting to solve the problem. We only do this as part of the presolver=automatic setting, if you tell the presolver that you want as much presolver as possible, you can use the presolver=aggressive setting. For your test cases here, that setting was still pretty fast and does the necessary reduction even in the 440x440 case. In general it is not advisable to use it since the presolver then might take several minutes even on reasonably small instances.

Rob's suggestion to implicitly substitute out the XIT variables is the way to go in your case.

I will see if I can come up with a way to make the presolver=automatic setting a little bit smarter in the future to get this improvement automatically.

If you encounter any other issue with the MILP solver, let us know.

Best regards

Philipp

ballardw · Posted 06-21-2021 01:45 PM

@genemroz wrote:

This is a follow-up question relate to my posting of yesterday. Rob Pratt responded with an excellent solution to my question that works well within the "basic" version that was posted with the question. But when I tried to implement his solution into a model version with more cameras, orientations and larger target areas, solution time for PROC OptModel becomes excessive. I didn't experience that problem with my earlier version of the optimization model. What's the best way to proceed?

Thanks,

Gene

Possibly "adjust expectations".

How "excessive" is the time? If you expect a complex model to run in 5 seconds and it takes 5 minutes I think that may just be optimism and not understanding just how many computations may be going on in the background and the additional elements are likely not adding computations in a simply linear manner.

Many "models" increase by orders of magnitude with each additional parameter/constraint imposed or exponentially.

If 3 cameras take, just as a number to discuss, 1,000,000,000 computations and adding one more multiplies that by 10, would you expect it to take the same amount of time? or roughly 10 times as long?

The other parameters such as "larger target size" could have a drastic impact on the number of calculations needed. Double the dimensions of an "area" and you square the size involved. The surface of a solid could be much greater than that.

SAS Innovate 2025: Register Now