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;
