I’ve been using the attached code (big hat tip to Rob Pratt) to optimize camera pointing in a network of meteor cameras. I’d like to extend this code to yield the second best and maybe even the top three solutions. I have no idea how to approach this. So I’m looking for guidance on how to make this happen.
I'm only including the code and no datasets because they're too large to readily attach. I'm hoping that the code and commenting is sufficiently clear that someone can take a look at it and get me pointed in the right direction.
I’ve occasionally had issues where my SAS ODA session is “Terminated Due to Excessive Resource Usage” when running this code. I definitely want to avoid that circumstance. So if termination is made significantly more likely by finding the top 2 or 3 solutions, then I’ll probably just stick with what is working.
Thanks in advance for any quidance you can provide,
Gene
/*SAS Procedure Optmodel is at the heart of the code. It analyzes the matrix of all possible orientations
for the cameras to be optimized together with the fields of view of the cameras with fixed
orientations. It seeks to find the orientations of the cameras to be optimized that
maximizes the number of gridpoints that are covered by 3 or more cameras */
Proc options group=(performance memory);
run;
Proc Optmodel printlevel=0;
options fullstimer;
/* Read CoverageMatrix data into index sets*/
set <str> stations;
read data work.stationdata4opt into stations=[station];
/*The set Station_Pairs statement below is used when two or more cameras are co-located.
It is used in conjuction with the set Azimuths statement below to forbid solutions
where adjacent azimuths are assigned to co-located cameras. Such a solution is
undesirable because because trajectory solutions from such pairs would have
a zero convergence angle. */
/* set Station_Pairs={<'US003G','US003S'>,<'US003G','US003T'>,<'US003S','US003T'>, */
/* <'US000L','US003G'>,<'US000L','US003S'>,<'US000L','US003T'>}; */
set <str> ORIENTATION_STRING=
/'45/35' '90/35' '135/35' '180/35' '225/35' '270/35' '315/35' '360/35'
'45/45' '90/45' '135/45' '180/45' '225/45' '270/45' '315/45' '360/45'
'45/55' '90/55' '135/55' '180/55' '225/55' '270/55' '315/55' '360/55'/;
set ORIENTATIONS = setof {s in orientation_string}
<input(scan(s,1,'/'),best.),input(scan(s,2,'/'),best.)>;
put orientations=;
set <str> Targets;
read data work.targets into TARGETS=[target];
num Matrix {stations, ORIENTATION_string, TARGETS};
read data work.RoO_Matrix into [station orientation_string target] Matrix=k;
/* Set Constants*/
/* Sets optimal number of cameras covering each target*/
%Let k=3;
num k=&k;
/*Declare Variables*/
var CHI {stations, ORIENTATIONS} binary;
var IsTriplyCovered {TARGETS} binary;
impvar XIT{target in TARGETS}=sum{station in stations, <a,e> in
orientations} matrix[station, a||'/'||e, target] *Chi[station,a,e];
/*Rob Pratt suggestion to improve execution time. For the
Fixed Cameras,all solutions are equivalent. This command
"fixes" the solution (CHI) for the Fixed Cameras ('USXXXX') at AZ=45 and ELV=35.
If there are no fixed cameras then code below should be commented out*/
fix CHI['USXXXX',45,35]=1;
/* It's also possible to "fix" the orientation of other cameras.
If it is known that a camera for optimization
has a blocked azimuth/elevation then set CHI = 0. See example below: */
/* fix CHI['US000L',90,35]=0; */
/* Declare Model*/
max NumTriplyCovered=sum{t in TARGETS} IsTriplyCovered[t];
/*Subject to Following Constraints*/
/*Each station can have only one optimal orientation */
con CHI_constraint {station in stations}:
sum{<a,e> in ORIENTATIONS}CHI[station,a,e] <=1;
/*Checks that target is triply covered...?*/
con TriplyCoveredCon {t in TARGETS}: XIT[t] >=k * IsTriplyCovered[t];
/*This constraint defines adjacent azimuths as azimuth_pairs. It then
checks that solution does not include any Station_Pair with adjacent azimuths*/
/* set AZIMUTHS = setof {<a,e> in ORIENTATIONS} a; */
/* put AZIMUTHS=; */
/* set AZIMUTH_PAIRS = {a1 in AZIMUTHS, a2 in AZIMUTHS: mod(a1-a2+360,360) <= 45}; */
/* put AZIMUTH_PAIRS=; */
/* */
/* con NotAdjacentAzimuths {<s1,s2> in STATION_PAIRS, <a1,a2> in AZIMUTH_PAIRS}: */
/* sum {<a,e> in ORIENTATIONS: a in {a1,a2}} (CHI[s1,a,e] + CHI[s2,a,e]) <= 1; */
/* expand; */
/* Call Solver and Saves Initial Results:*/
/*relobjgap (0-1) specifies stopping criterion based on the best
integer objective and the best bound on the objective function value for
the Initial Result.
After some experimentation with a 5 million observation dataset,
I've specified 0.9 as having the best performance when used in conjunction
with subsequent optimization code below. With smaller datasets,
a smaller value could be used. Some experimentation may be required.*/
solve with milp / symmetry=moderate presolver=aggressive relobjgap=.9;
create data work.InitialSolution_RoO from
[station a e]={c in stations, <a,e> in orientations: CHI[c,a,e]>0.5} CHI;
/* Takes Initial Solution as starting point for additional optimization on a pair by pair basis.
This became important when optimizing many (>4-5?) stations because cpu times were becoming excessive.
This approach, although a bit of a short-cut, is computationally more efficient.*/
num bestObj;
bestObj = _OBJ_;
num improved;
fix Chi;
num numSubsets init 0;
set SUBSETS = 1..numSubsets;
set <str> STATIONS_s {SUBSETS};
for {c1 in STATIONS, c2 in STATIONS diff {c1}: c1 < c2} do;
numSubsets = numSubsets + 1;
STATIONS_s[numSubsets] = {c1,c2};
end;
do until (improved = 0);
improved = 0;
for {s in SUBSETS} do;
put STATIONS_s[s]=;
for {c in STATIONS_s[s], <a,e> in orientations} unfix Chi[c,a,e];
solve with milp / primalin;
if bestObj < _OBJ_ then do;
bestObj = _OBJ_;
improved = 1;
end;
fix Chi;
end;
end;
unfix Chi;
create data work.OptimalSolution_RoO from
[station a e]={c in stations, <a,e> in orientations: CHI[c,a,e]>0.5} CHI;
quit;
/*Print Results*/
proc print data=work.initialsolution_RoO;
Title "Initial Solution";
proc print data=work.optimalsolution_RoO;
Title "Optimal Solution";
run;
/* Create Macro Variables for Station and Optimal Solution*/
data macrovariables;
set optimalsolution_RoO;
suffix=put(_n_,5.);
call symput (cats('Sta_Code',suffix),Station);
call symput (cats('Opt_Orientation',suffix),catx('/',a,e));
run;
Thanks in advance for any guidance you can provide.
... View more
