BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
genemroz
Quartz | Level 8

All,

I've been working to develop an optimization model to optimize the aiming of multiple cameras to cover a three-dimensional Cartesian target area within a three-dimensional "regional" grid. The objective function of the algorithm I've implemented optimizes among 8 possible azimuth angles and 3 elevation angles for each of up to 24 cameras. The objective function in Proc OPTMODEL maximizes the count of grid points covered by the cameras within the target area. Furthermore, while each gridpoint can be covered by multiple cameras, more than three cameras per gridpoint do not inrease the objective value. The current optimization model gives the same weight to bringing a gridpoint from zero coverage to single-camera coverage as it does for bringing a gridpoint from single- to double- or double- to triple-camera coverage.

 

However, for my purposes, bringing a gridpoint from zero coverage to single coverage is of no value. I seek to optimize the overlap of the camera fields of view by maximizing the number of gridpoints that are triply and, secondarily, doubly covered.

 

I would like to modify the objective function so that it maximizes the number of gridpoints covered by three cameras. As above, it is acceptable if gridpoints are covered by more than three cameras, but there should be no additional increase in the objective value for this condition.

 

I'd be very grateful is someone can provide guidance on how to change the optimization function or constraints to achieve my goal.

 

The code below is a "basic" version of the optimization model but with three cameras fixed at a single elevation angle and with 8 available azimuth possibilities. I've also limited the target area to a single altitude level. If you run the code and examine the plot of camera coverage, I think it's pretty obvious that a better solution would be for camera US001E to point northeast. So that would be the desired outcome of the revised optimization model.

 

Thanks in advance for any guidance you can provide.

 

Gene

 

/*Basic Camera Network Optimization Model*/
data 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
US001E 33.656305 -117.560879 .439 135 45
;
run;
data stationdataexisting;
	set stationdataexisting;
	/*Define station grid coordinates relative to Needles CA*/
	xpos=(long-(-114.623))*(111.321*cos(constant('pi')/180*lat));
	ypos=(lat-(34.766))*111.19;
	zpos=alt;
run;

/*Define 1 single elevation angle for each camera in radians. If a certain elevation
is to be omitted, enter a zero to indicate missing*/
options mprint symbolgen;
%LET LO=34;
/* Enter dimension of one side of Region of Optimization Matrix (RoO) centered on Needles CA*/
%let RoO=500;

/*Define camera orientation unit vectors*/
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;
	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);

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

	/*EAST AZIMUTH*/
	/*AZEA_ALTLO*/
	ealo_xhat=sin(azea);
	ealo_yhat=cos(azea);

	/*SOUTHEAST AZIMUTH*/
	/*AZSE_ALTLO*/
	selo_xhat=sin(azse);
	selo_yhat=cos(azse);

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

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

	/*WEST AZIMUTH*/
	/*AZWE_ALTLO*/
	welo_xhat=sin(azwe);
	welo_yhat=cos(azwe);

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

	/* Create Coverage Matrix*/
data work.CoverageMatrix;
	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(24) _temporary_;
	set work.orientationvectors;

	/*North Unit Vectors*/
	vector(1)=nlo_xhat;
	vector(2)=nlo_yhat;
	vector(3)=zhatlo;

	/*Northeast Unit Vectors*/
	vector(4)=nelo_xhat;
	vector(5)=nelo_yhat;
	vector(6)=zhatlo;

	/*East Unit Vectors*/
	vector(7)=ealo_xhat;
	vector(8)=ealo_yhat;
	vector(9)=zhatlo;

	/*Southeast Unit Vectors*/
	vector(10)=selo_xhat;
	vector(11)=selo_yhat;
	vector(12)=zhatlo;

	/*South Unit Vectors*/
	vector(13)=slo_xhat;
	vector(14)=slo_yhat;
	vector(15)=zhatlo;

	/*Southwest Unit Vectors*/
	vector(16)=swlo_xhat;
	vector(17)=swlo_yhat;
	vector(18)=zhatlo;

	/*West Unit Vectors*/
	vector(19)=welo_xhat;
	vector(20)=welo_yhat;
	vector(21)=zhatlo;

	/*Northwest Unit Vectors*/
	vector(22)=nwlo_xhat;
	vector(23)=nwlo_yhat;
	vector(24)=zhatlo;

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

	do Ty=Tymin to Tymax by Interval;

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

			/* 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 8  possible orientations of station c
				to determine whether Target(Tz, Ty, Tx) is covered (k=1) or not (k=0)*/
				do j=1 to 22 by 3;

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

					if j=4 then
						orientation='NELO';

					if j=7 then
						orientation='EALO';

					if j=10 then
						orientation='SELO';

					if j=13 then
						orientation='SOLO';

					if j=16 then
						orientation='SWLO';

					if j=19 then
						orientation='WELO';

					if j=22 then
						orientation='NWLO';

					/*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'));
					UpperVertFoVlimit=arsin(zhat)*(180/constant('pi'))+22;
					LowerVertFoVlimit=arsin(zhat)*(180/constant('pi'))-22;

					/*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'));

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

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

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

							/*end 1*/
						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;
	stop;
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;
	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=/'NOLO' 'NELO' 'EALO' 'SELO' 'SOLO' 'SWLO' 'WELO' 
		'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*/
	/*number of stations*/
	%Let N=3;
	num N=&N;

	/* optimal number of cameras covering each target*/	
	%Let k=3;
	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;

	/*Sadik Equation 9. I set k=3 to signify maximize up to 3 cameras per target.*/
	var PSI{TARGETS} integer <=k;

	/* Declare Model*/
	/*Sadik Equation 7--the function to be maximized*/
	max OptCoverage=sum{target in TARGETS} Psi[target];

	/*Subject to Following Constraints*/
	con CHI_constraint {station in stations}: 
	sum{orientation in ORIENTATIONS}CHI[station, orientation] <=1;

	/* Sadik Equation 8*/
	/* Sadik: Syntax for expression XIT and  XIT_N for use in Sadik Equation 8 below*/
	var XIT{TARGETS} integer;
	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 from 
	[station orientation]={c in stations, o in orientations: CHI[c, o]>0.5} CHI;
quit;

data work.optimalsolution;
	set work.optimalsolution;

	if orientation='NOLO' then
		do Az=0;
			Alt=&lo;
		end;
	else
		;

	if orientation='NELO' then
		do Az=45;
			Alt=&lo;
		end;
	else
		;

	if orientation='EALO' then
		do Az=90;
			Alt=&lo;
		end;
	else
		;

	if orientation='SELO' then
		do Az=135;
			Alt=&lo;
		end;
	else
		;

	if orientation='SOLO' then
		do Az=180;
			Alt=&lo;
		end;
	else
		;

	if orientation='SWLO' then
		do Az=225;
			Alt=&lo;
		end;
	else
		;

	if orientation='WELO' then
		do Az=270;
			Alt=&lo;
		end;
	else
		;

	if orientation='NWLO' then
		do Az=315;
			Alt=&lo;
		end;
run;

proc sort data=work.OptimalSolution;
	by station;
run;

/*Create macro to select Region of Optimization (RoO) dimensions*/
/* Enter dimension(km)of one side of square centered on ABQ below:*/
options mprint symbolgen;

%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;
	by station;
run;

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

Data work.stationdataoptimal;
	merge work.stationdataexisting work.optimalsolution;
	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;
	Interval=10;
	Target=0;
	Tz=100;

	/* 	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);

				/*Calculate Elevation angle and  vertical FOV limits between Camera Location and Target Location*/
				Elev2Target=arsin(CTz/CTdist)*(180/constant('pi'));
				UpperVertFoVlimit=arsin(xyz(5, c))*(180/constant('pi'))+22;
				LowerVertFoVlimit=arsin(xyz(5, c))*(180/constant('pi'))-22;

				/*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*/
				/*1*/
				if 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 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;
	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;

data citycoord;
	input city & $15. lat long;

	/*Be sure to follow city name with TWO SPACES before lat entry*/
	datalines;
+ LAX  33.942 -118.408
+ Las Vegas  36.080 -115.152
+ San Diego  32.734 -117.19
+ St. George  37.036 -113.51
+ Needles  34.766 -114.623
+ Yuma  32.669 -114.599
+ Bakersfield  35.3733 -119.0187
+ Victorville  34.5362 -117.2928
+ Barstow  34.8958 -117.0173
+ Death Valley  36.5323 -116.9325
;
run;

data city_anno_Needles (drop=city lat long xpos ypos);
	set citycoord;
	length function $8;

	/*location of Needles*/
	xpos=(long-(-114.623))*91.09;
	xpos=round(xpos, 10);
	ypos=(lat-(34.766))*111.19;
	ypos=round(ypos, 10);
	label=city;
	x1=xpos;
	y1=ypos;
	anchor='left';
	x1space='datavalue';
	y1space='datavalue';
	textsize=8;
	widthunit='data';
	width=200;
	function='text';
run;

data stationcoord;
	input station & $15. lat long;

	/*Be sure to follow station name with TWO SPACES before lat entry*/
	datalines;
+ US000S  34.291 -116.384
+ US000V  33.960 -117.259
+ US001E  33.656 -117.561
;
run;

data station_Needles_anno;
	set stationcoord;
	length function $8;

	/*Define Center of Grid as location of Needles*/
	xpos=(long-(-114.623))*91.09;
	xpos=round(xpos, 10);
	ypos=(lat-(34.766))*111.19;
	ypos=round(ypos, 10);
	label=station;
	x1=xpos;
	y1=ypos;
	anchor='left';
	x1space='datavalue';
	y1space='datavalue';
	textsize=8;
	widthunit='data';
	width=200;
	function='text';
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, station, poly and text annotations*/
data anno_all;
	set anno_poly anno_text city_anno_needles station_Needles_anno;
run;

/*Build of HEATMAPPARM Code*/
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;

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 "degrees";
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;

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";
	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;
1 ACCEPTED SOLUTION

Accepted Solutions
RobPratt
SAS Super FREQ

To maximize the number of targets that are covered by at least three cameras, replace Psi, the original objective, and the XIT_CON2 and XIT_CON3 constraints with the following:

   /* maximize number of triply covered targets */
   var IsTriplyCovered {TARGETS} binary;
   max NumTriplyCovered = sum {t in TARGETS} IsTriplyCovered[t];
   con TriplyCoveredCon {t in TARGETS}:
      XIT[t] >= 3 * IsTriplyCovered[t];

In SAS Viya, you can alternatively use an indicator constraint as follows:

   con TriplyCoveredCon {t in TARGETS}:
      IsTriplyCovered[t] = 1 implies XIT[t] >= 3;

View solution in original post

4 REPLIES 4
RobPratt
SAS Super FREQ

To maximize the number of targets that are covered by at least three cameras, replace Psi, the original objective, and the XIT_CON2 and XIT_CON3 constraints with the following:

   /* maximize number of triply covered targets */
   var IsTriplyCovered {TARGETS} binary;
   max NumTriplyCovered = sum {t in TARGETS} IsTriplyCovered[t];
   con TriplyCoveredCon {t in TARGETS}:
      XIT[t] >= 3 * IsTriplyCovered[t];

In SAS Viya, you can alternatively use an indicator constraint as follows:

   con TriplyCoveredCon {t in TARGETS}:
      IsTriplyCovered[t] = 1 implies XIT[t] >= 3;
genemroz
Quartz | Level 8

Rob,

Thank so much for your quick response and suggestion.

 

I implemented your suggestion (as I understood it) but get an error message that I don't understand.  The error message says symbol 'XIT' is unknown.   The log of the offending code is below.   Perhaps I didn't incorporate your suggestion correctly.  Thanks again for your help with this

 

Gene

368        Proc optmodel;
 369        /* Read CoverageMatrix data into index sets*/
 370        set <str> stations;
 371        read data work.stationdataexisting into stations=[station];
 NOTE: There were 3 observations read from the data set WORK.STATIONDATAEXISTING.
 372        set <str> ORIENTATIONS=/'NOLO' 'NELO' 'EALO' 'SELO' 'SOLO' 'SWLO' 'WELO'
 373        'NWLO'/;
 374        set <num> TARGETS;
 375        read data work.targets into TARGETS=[target];
 NOTE: There were 2601 observations read from the data set WORK.TARGETS.
 376        num Matrix {stations, ORIENTATIONS, TARGETS};
 377        read data work.CoverageMatrixOpt into [station orientation target] Matrix=k;
 NOTE: There were 62424 observations read from the data set WORK.COVERAGEMATRIXOPT.
 378        
 379        /* Set Constants*/
 380        /*number of stations*/
 381        %Let N=3;
 382        num N=&N;
 SYMBOLGEN:  Macro variable N resolves to 3
 383        
 384        /* optimal number of cameras covering each target*/
 385        %Let k=3;
 386        num k=&k;
 SYMBOLGEN:  Macro variable K resolves to 3
 387        
 388        /*Declare Variables*/
 389        /* Sadik Equation 11--Sets CHI as binary variable that is 1 when using station i and POSITION j else 0.*/
 390        var CHI {stationS, ORIENTATIONS} binary;
 391        
 392        /*Sadik Equation 9. I set k=3 to signify maximize up to 3 cameras per target.*/
 393        /* var PSI{TARGETS} integer <=k; */
 394        
 395        
 396        /* Declare Model*/
 397        /*Sadik Equation 7--the function to be maximized*/
 398        /* max OptCoverage=sum{target in TARGETS} Psi[target]; */
 399        
 400        /*Rob Pratt's suggested function to be maximized*/
 401        var IsTriplyCovered {TARGETS} binary;
 402        max NumTriplyCovered = sum{t in TARGETS} IsTriplyCovered[t];
 403        con TriplyCoveredCon {t in TARGETS}: XIT[t] >=3 * IsTriplyCovered[t];
                                                  ___
                                                  525
 ERROR 525-782: The symbol 'XIT' is unknown.
 
 404        
 405        /*Subject to Following Constraints*/
 406        con CHI_constraint {station in stations}:
 407        sum{orientation in ORIENTATIONS}CHI[station, orientation] <=1;
 408        
 409        /* Sadik Equation 8*/
 410        /* Sadik: Syntax for expression XIT and  XIT_N for use in Sadik Equation 8 below*/
 411        var XIT{TARGETS} integer;
 412        con XIT_CON1{target in TARGETS}:
 413           XIT[target]=sum{station in stations, orientation in ORIENTATIONS}
 414        matrix[station, orientation, target] * Chi[station, orientation];
 415        /* Removed per Rob Pratt's suggested modification*/
 416        /* con XIT_CON2{target in TARGETS}:XIT[target]/N <=PSI[target]; */
 417        /* con XIT_CON3{target in TARGETS}:PSI[target] <=XIT[target]; */
 418        
 419        /* Call Solver and Save Results*/
 420        solve;
 NOTE: Problem generation will use 2 threads.
 NOTE: Previous errors might cause the problem to be resolved incorrectly.
 ERROR: The constraint 'TriplyCoveredCon' has an incomplete declaration.
 NOTE: The problem has 5226 variables (2601 free, 0 fixed).
 NOTE: The problem has 2625 binary and 2601 integer variables.
 NOTE: The problem has 2604 linear constraints (3 LE, 2601 EQ, 0 GE, 0 range).
 NOTE: The problem has 10042 linear constraint coefficients.
 NOTE: The problem has 0 nonlinear constraints (0 LE, 0 EQ, 0 GE, 0 range).
 NOTE: The OPTMODEL presolver is disabled for linear problems.
 NOTE: Unable to create problem instance due to previous errors.
 421        create data work.OptimalSolution from
 422        [station orientation]={c in stations, o in orientations: CHI[c, o]>0.5} CHI;
 NOTE: The data set WORK.OPTIMALSOLUTION has 0 observations and 3 variables.
 423        quit;
 NOTE: The SAS System stopped processing this step because of errors.
 NOTE: PROCEDURE OPTMODEL used (Total process time):
       real time           0.16 seconds
       user cpu time       0.17 seconds
       system cpu time     0.02 seconds
       memory              15612.78k
       OS Memory           51272.00k
       Timestamp           06/20/2021 06:27:22 PM
       Step Count                        204  Switch Count  6
       Page Faults                       0
       Page Reclaims                     4351
       Page Swaps                        0
       Voluntary Context Switches        70
       Involuntary Context Switches      1
       Block Input Operations            0
       Block Output Operations           280
       
 424        
RobPratt
SAS Super FREQ
The new constraint involves XIT and so must appear after that variable is declared.
genemroz
Quartz | Level 8

Of course!  Your suggestion made a huge improvement to the results.   I'm marking it as Accepted.

 

Thanks again!

 

Gene

SAS Innovate 2025: Register Now

Registration is now open for SAS Innovate 2025 , our biggest and most exciting global event of the year! Join us in Orlando, FL, May 6-9.
Sign up by Dec. 31 to get the 2024 rate of just $495.
Register now!

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.

Discussion stats
  • 4 replies
  • 764 views
  • 2 likes
  • 2 in conversation