Approach for developing a warm start for Proc Optmodel?

[genemroz]

Esteemed Advisers:

I'm looking for some guidance on the best approach(es) to developing a Warm Start solution for use in conjunction with Proc Optmodel and the PRIMALIN option. Is it appropriate to use, say, Proc LP? What approach would you recommend?

Thanks,

Gene

[RobPratt]

Is this for the MILP solver?  Are you able to share your code and data?

A warm start solution can come from anywhere, but I would not suggest using PROC LP, which is a legacy procedure that is no longer actively maintained and was last documented in SAS/OR 14.1 (released in 2015).

Often, you can construct a good warm start solution by solving a simplified version of your original problem, and for that a good choice is to use PROC OPTMODEL itself, with multiple SOLVE statements in the same OPTMODEL session.

[genemroz]

Rob,

Thanks for the prompt reply.  I'll try your suggestion and will prepare code and data files for posting.

Thanks,

Gene

[genemroz]

Rob:

Thanks for offering to help.

The basic issue here is that this optimization problem is taking so long to solve that I keep getting unceremoniously kicked off of the SAS OnDemand for Academics server (a warning would be nice) for "excessive resource usage".

You helped me late last year with a similar problem by suggesting the invocation of the SYMMETRY option in Proc Optmodel. It helped speed the solution in that particular instance of the problem but, in this different instance, invoking the SYMMETRY option doesn't result in any orbits being found and appears to have no benefit.

Recently I came across a paper by Klotz and Newman entitled "Practical Guidelines for Solving Difficult Mixed Integer Linear Programs"--just what I need! In that paper they note that "lack of progress in the best integer solution" (which is what I seem to be experiencing) is best addressed by providing an initial solution. So that's what got me interested in using the PRIMALIN option and how I came to ask about the best approach to finding a Warm Start solution.

It puzzles me why this problem is so difficult to solve. In my uneducated opinion it seems to be a relatively straightforward (set cover or assignment?) problem: find the best orientation for 9 cameras that maximizes three-camera coverage over a set of Cartesian gridpoints. Each camera has only 16 possible orientations so there are only 9x16=144 possible solutions to choose from. How hard can that be?

I've prepared and attached the following the files:

1 Dataset (FoV.zip) that contains the field of view for each of 16 possible camera orientations. This is a zipped .csv file that will have to be un-zipped and imported into a SAS dataset.

2. Code (Data Preparation.sas) that invokes dataset (FoV) to prepare the "Region of Optimization Matrix" (RoO_Matrix) dataset. This is a 3 dimensional matrix of 9 stations (cameras), 16 orientations and several thousand targets (depending on the size of the Region of Optimization). Each element of the matrix is binary where k=1 when a target is covered by (i.e. within the field of view of) a given combination of camera and orientation or, otherwise, k=0.

Note: the path to dataset FoV needs to be specified in this Code at line 59.

3. The Optimization code (ProcOptodel.sas) consists of three parts:

Part 1 invokes Proc Optmodel, creates index sets for stations, orientation_string and targets, reads dataset RoO_Matrix, defines the model and solves it using MILP to provide the Warm Start solution by fixing the orientations of 4 stations based on my best guess. With RELOBJGAP=0.0, this problem solves in about 90 seconds.

Part 2 invokes Proc Optmodel using MILP with the PRIMALIN option to optimize all 9 stations. With RELOBJGAP=0.25 or lower the solution time exceeds 10 minutes, Best Integer progress is slows to a crawl and is at risk of getting a "Session terminated due to excessive resource usage" message when run on SAS OnDemand for Academics server.

Part 3. Produces a report and graphic output.

I would like for the code to solve without having to accept the sub-optimal solution that results when setting RELOBJGAP to a value that will stop the execution before the session is "terminated due to excessive resource usage".

Let me know if you have any questions or problems accessing the attached files.

Thanks so much for offering to help,

Gene

[RobPratt]

I'll take a look soon, but first a couple of quick comments regarding "Each camera has only 16 possible orientations so there are only 9x16=144 possible solutions to choose from. How hard can that be?":

1. If the orientation for each camera can be chosen independently, there are 9x16 decision variables but 16^9 solutions.

2. The number of possible solutions is not necessarily a good indicator of difficulty.  An (n+1)-city traveling salesman problem has n! solutions (or n!/2 if undirected) and is difficult for large n.  On the other hand, an nxn linear assignment problem also has n! solutions but is easy even for large n.

[RobPratt]

I am getting a WARNING from READ DATA:

1256 read data work.targets into TARGETS=[target];
NOTE: There were 12617 observations read from the data set WORK.TARGETS.
1257 num Matrix {stations, ORIENTATION_string, TARGETS};
1258 read data work.RoO_Matrix into [station orientation_string target]
1259 Matrix=k;
WARNING: 50832 non-missing values were discarded due to invalid target indices.
NOTE: There were 1629694 observations read from the data set WORK.ROO_MATRIX.

There are 12970 distinct values of target in WORK.ROO_MATRIX but only 12617 observations in WORK.TARGETS.

This discrepancy then yields these ERRORs:

1307 solve with milp/logfreq=60 relobjgap=0.0;
NOTE: Problem generation will use 4 threads.
ERROR: The array element 'Matrix[WELLIN,'45/35','90/90/12']' has no value at line 1270 column 36.
ERROR: The array element 'Matrix[WELLIN,'45/35','90/90/11']' has no value at line 1270 column 36.
ERROR: The array element 'Matrix[WELLIN,'45/35','90/90/10']' has no value at line 1270 column 36.
ERROR: The array element 'Matrix[WELLIN,'45/35','90/60/12']' has no value at line 1270 column 36.
NOTE: Unable to create problem instance due to previous errors.

Please double check the data.  It might also help if you can provide the exact code you used to convert from csv to SAS data set.

[genemroz]

Rob,

I apologize for this situation. I'm trying to figure out what went wrong. I can see from the ERROR log that the last element in the array '90/90/12' should read '90'90'120', so somehow the last digit got truncated.  The csv file I attached was created by exporting the SAS dataset and then I zipped it.  Maybe that wasn't a good idea?  I've attached an uncompressed version of FoV.csv.  Is there way to send the SAS dataset directly?

My apologies,

Gene

[RobPratt]

The csv itself seems to be OK, and after clearing WORK and trying again, I no longer see the error.  Just to confirm, did you get objective value 1102 for the first solve?

[genemroz]

Yes. Objective Value of 1102 for the first solve with 4 stations with fixed orientations.

Gene

[RobPratt · Accepted solution]

Here's a "MILP local search" approach to improve the 1102 solution to 1207 quickly by calling the solver in a loop.  The idea is to fix the orientations for all but a small subset of the stations and optimize over the resulting neighborhood of the currently best integer feasible solution.  Subsets of size 2 performed the best here, but I left the code for subsets of size 1 or 3 commented out for reference.  The DO UNTIL loop ends when the current solution is locally optimal in the sense that changing only two stations cannot improve the solution.

num bestObj;
bestObj = _OBJ_;
num improved;
fix Chi;
num numSubsets init 0;
set SUBSETS = 1..numSubsets;
set <str> STATIONS_s {SUBSETS};
/*for {c in STATIONS} do;*/
/*   numSubsets = numSubsets + 1;*/
/*   STATIONS_s[numSubsets] = {c};*/
/*end;*/
for {c1 in STATIONS, c2 in STATIONS diff {c1}: c1 < c2} do;
   numSubsets = numSubsets + 1;
   STATIONS_s[numSubsets] = {c1,c2};
end;
/*for {c1 in STATIONS, c2 in STATIONS diff {c1}, c3 in STATIONS diff {c1,c2}: c1 < c2 < c3} do;*/
/*   numSubsets = numSubsets + 1;*/
/*   STATIONS_s[numSubsets] = {c1,c2,c3};*/
/*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;

If you want, you can then call the solver again with PRIMALIN to solve the original problem with no fixed variables.  After running for a long time, the solver did not improve BestInteger but improved BestBound to 1254, so 1207 is within 4% of optimal and might actually be optimal.

[genemroz]

Rob,

Thank you so much for taking the time to help with this problem.

I've incorporated your 'MILP local search" into my code and it ran quickly and successfully.  I will mark your solution as accepted.

I find it curious but reassuring that I got the same objective value as you (1207) when I used the 'solve with lp relaxint' statement for the first solve, followed by a 'solve with milp/primalin' .  However, this approach required nearly 20 minutes to achieve that objective value with a still large gap (20%) and I wasn't very confident of it's optimality.

I will implement your approach for the several other instances of this problem that I have on my plate to address.  I hope it will prove robust enough to handle those as well.  Are there any caveats that I should be aware of when applied to other problems of this nature?

Thanks again for your help with this.   This project would not be able to go forward without your assistance.

Gratefully,

Gene

[RobPratt]

I'm always glad to help, and I appreciate your questions in this forum, where other users can benefit from the responses.

In reply to your request for caveats, I should just mention that my suggestion is an improvement heuristic and so there is no guarantee of optimality or even improvement.  But, as you have seen, this "MILP local search" idea can sometimes be useful for exploiting problem structure to find better feasible solutions quickly.

[genemroz]

Yes, Objective Value of 1102 for the first solve.

Gene