Showing posts with label AMPL. Show all posts
Showing posts with label AMPL. Show all posts

Wednesday, May 23, 2012

K Best Solutions in AMPL/CPLEX

This is a continuation of my earlier post on how to find the K best solutions to a mixed integer linear program. Per a reader request, I'll show how to do it using AMPL+CPLEX. Note that cplexamp is the name of the CPLEX executable that interfaces with AMPL.  The example is that same one from the previous post, a binary knapsack model.

Here's the AMPL code:

#
# K best solutions using AMPL + CPLEX
# Test problem is a binary knapsack
#
# model parameters
#
param n default 20;  # problem size
param k default 3;   # solutions to find
param weight {i in 1..n} := i*sqrt(i);  # item weights
param value {i in 1..n} := i^2;  # item values
param capacity := (2/3)*sum {i in 1..n} weight[i];  # knapsack capacity
#
# parameters controlling the solution process
#
param method in {1,2} default 1;  # method:
                                  # 1=solve MIP, then populate;
                                  # 2=populate only
param intensity in 0..4 default 3;  # CPLEX pool intensity
param limit default 20;     # number of pool solutions to generate
#
# the binary knapsack model
#
var x {1..n} binary;
maximize TotalValue: sum {i in 1..n} value[i]*x[i];
s.t. TotalCapacity: sum {i in 1..n} weight[i]*x[i] <= capacity;
#
# print out solution parameters
#
print "Problem size = ", n, ", seeking ", k, " solutions, method = ",
                      method, ", intensity = ", intensity, ", limit = ",
                      limit;
#
# set up CPLEX
#
option solver cplexamp;
option cplex_options ("populate=" & method & 
                      " poolstub=kbpool poolcapacity=" & k & 
                      " poolintensity=" & intensity & 
                      " poolreplace=1 populatelim=" & limit);
#
# solve the model and show the pool solutions
#
solve;
for {i in 1 .. Current.npool} {
  solution ("kbpool" & i & ".sol");
  display _varname, _var;
}

Most of the CPLEX options can be inferred from the AMPL parameters assigned to them.  There are two I should probably explain:
  • poolstub is a file name that CPLEX will use to store each of the solutions in the pool (so in my case I end up with files kbpool1.sol, kbpool2.sol, and kbpool3.sol), which the final loop reads back in; and
  • poolreplace tells CPLEX how to decide which solution to boot out of the pool when the pool is full and a new solution is found (value 1 means keep the solutions with the best objective value; other choices are to keep the most recently found solutions or to strive for diversity).
One last note: the .npool suffix (here attached to the current problem) returns the number of solutions in the solution pool at termination (which equals the number of .sol files generated).

Addendum: The last note is apparently not the last note. :-) Guido's comment below reminded me that an important point from the previous post needs to be repeated here (particularly if someone stumbles on this without reading the prior post).  The title of the post notwithstanding, the solution pool approach may not actually produce the K best solutions, unless the population limit is set quite high and the pool intensity is at one of the higher settings.  Both those settings changes are likely to increase solution times.  As was seen in the previous post on this, the moderate settings used in the code example produce several "good" solutions, one of which is optimal, but they do not always produce the second best, third best etc. solutions.

Addendum #2: To answer a question below, we can store the various solutions in a parameter (xx) in the script by adding the following code to the end of the script:

#
# store the solutions
#
param xx {i in 1 .. Current.npool, j in 1 .. n};
for {i in 1 .. Current.npool} {
  for {j in 1 .. n} {
    solution ("kbpool" & i & ".sol");
    let xx[i, j] := _var[j];
  }
}
display xx;  # verify what was stored

Sunday, May 29, 2011

Peeking Inside the Black Box

From the perspective of an analyst, rather than a developer of algorithms, the computational side of operations research might be described in terms of a black box: insert model and data; turn crank; make sense of the results. When it comes to the question of how much the analyst needs to know about what goes on inside the box, I find myself walking a somewhat narrow line between the Um and the Yang.

When teaching statistics to end users, I resist the exhortations of some of my colleagues that students must do computations by hand to “really learn” about statistics. For instance, said colleagues will argue that in order for a student to grasp the concept of correlation, they must take a set of observations of a pair of variables and (manually) center the variables, compute squares and cross-products, total the squares and cross-products, and use the totals to compute standard deviations and the Pearson correlation coefficient. To me, that is just gratuitous grunt work. Students need to know the domain of the correlation coefficient (-1 to +1), when to use it (looking for linear relationships) or not (looking for nonlinear relationships), and its interpretation (including the dangers of making a leap from correlation to causation). The actual computation can safely be left to software.

On the other hand, the widespread availability of modeling software and powerful solvers has changed the practice of optimization in a very fundamental way. Once the exclusive province of people who had taken extensive coursework in the mathematical theory behind the pertinent algorithms (not to mention formal training in numerical analysis), optimization is now increasingly being undertaken by people who apparently have minimal formal training, and who may understand the contextual domain of the problem they are trying to solve but perhaps do not understand the mathematical underpinnings of the methods being used. I've seen evidence of this in a number of forums recently. One of the symptoms is a failure to understand issues relating to rounding error and numerical stability, ranging from the subtle (poor scaling leading to incorrect answers or solver failures) to the blatant (“My binary variable came out 0.99999999999 -- what went wrong??”).

A few recent posts suggested that some users either do not understand the difference between local and global optima or do not understand the role of convexity in optimization models. In particular, I cringe whenever I see nonlinear equality constraints in a model, whereas they don't seem to bother the model authors very much. (I'm not saying that nonlinear equality constraints are fundamentally wrong; I'm just saying that I'd really like not to have any in my model.) A very simple example will illustrate why.

Consider the following trivial optimization problems:\[ \begin{array}{clccclc} \textrm{(D)} & \min & x & \,\,\, & \textrm{(C)} & \min & x\\ & \textrm{s.t.} & x^{2}+y^{2}\le1 & & & \textrm{s.t.} & x^{2}+y^{2}=1\\ & & -2\le x\le2 & & & & -2\le x\le2\\ & & -2\le y\le2 & & & & -2\le y\le2\end{array}\]The (convex) feasible region of problem (D) is a disk, while the (nonconvex) feasible region of problem (C) is a circle. The solution to both problems is clearly $(x,y)=(-1,0)$.

Let's perform a little experiment (using AMPL to model the problems and MINOS 5.51 to solve them). We will solve both (D) and (C) with initial values $(x,y)=(cos(\theta),sin(\theta))$ where $\theta$ is an integer multiple of $\frac{\pi}{4}$. The AMPL code is as follows:

# problem definitions
var x >= -2 <= 2;
var y >= -2 <= 2;
minimize Obj: x;
s.t. Disk: x^2 + y^2 <= 1;
s.t. Circle: x^2 + y^2 == 1;
problem D: x, y, Obj, Disk;
problem C: x, y, Obj, Circle;

# bit used to define starting locations
param pi4 := atan(1); # pi/4
set ANGLE := 0..7; # from 0 to 7/4*pi

# storage for optimal solutions
param dsol {ANGLE};
param csol {ANGLE};

# solve disk problem from all starting angles
problem D;
for {t in ANGLE} {
  let x := cos(t*pi4);
  let y := sin(t*pi4);
  solve;
  let dsol[t] := x;
}

# solve circle problem from all starting angles
problem C;
for {t in ANGLE} {
  let x := cos(t*pi4);
  let y := sin(t*pi4);
  solve;
  let csol[t] := x;
}

display dsol, csol;
The results? For all starting angles, MINOS determines that $(x,y)=(-1,0)$ is the optimal solution to (D). That also applies for (C), with one glaring exception: when $\theta=0$, MINOS decides that the initial value $(x,y)=(1,0)$ is the optimal solution.

The figure below tells the tale. MINOS uses a reduced gradient method, which is a souped-up form of steepest descent. If allowed to, it will move in the direction pointed to by the gradient of the objective function at the current point when maximizing, and in the opposite direction when minimizing. (In this example, the gradient is a constant $(1,0)$.) When movement in that direction is infeasible, it will project the positive or negative gradient onto the hyperplane tangent to the feasible region at the current point, and search in that direction.


The red arrows indicate the (constant) negative gradient direction; the blue arrows are the projections of the gradient onto the unit circle. Having a hard time seeing the blue arrows anchored at $(\pm1,0)$? That's because they have length zero. So for any starting point other than $(\pm1,0)$ in problem (D), MINOS moves right to left across the unit disk until it reaches the boundary (the unit circle), then wends its way around the circle until it reaches the optimum at  $(-1,0)$. If it starts from $(+1,0)$, it cuts straight across the $x$-axis and arrives directly at the optimum. For problem (C), cutting across the disk is not an option. Given any start other than   $(\pm1,0)$, MINOS winds around the circle. Started at $(+1,0)$, however, MINOS sees a zero projection of the gradient, declares victory, stages a photo op and goes home.

The behavior of MINOS on problem (C) is correct given what's inside the black box. MINOS only guarantees local optima, and even that can only be delivered if the problem is moderately well behaved and/or the starting point is chosen well. Ipopt, another popular nonlinear solver, has similar limitations. From some recent posts that I've seen, though, I have the impression that some users are creating highly nonconvex models and expecting optimal solutions they are not destined to get. Compounding this, of course, is the fact that the solver will generally not know it has arrived at a suboptimal solution, and therefore will not alert the user that the solution is suboptimal.

Perhaps modeling languages and solvers should ship with a “warning label” stating that use of the software without first reading a book or taking a course in linear/nonlinear/integer programming is inadvisable.

Speaking of books, I can recommend David Luenberger's "Linear and Nonlinear Programming" (an updated version of which has apparently been published with coauthor Yinyu Ye).

Sunday, March 20, 2011

Semicontinuous Variables

This question pops up from time to time on various forums.  In a math programming model, how does one handle a variable that can either be zero or in some range not containing zero?  The mathematical description of the variable would be something like $$x\in {0} \cup [a,b]$$where $a>0$.  A couple of common applications in supply chain models include quantity discounts (a vendor offers a discounted price if you buy in bulk; $x$ is the amount purchased at the discounted price and $a$ is the minimum amount to qualify for the discount) and minimum lot sizes (where $x$ is the amount of a custom-made product obtained from a vendor who requires a minimum order of $a$ to justify the setup cost of the order).

The variable $x$ is known in math programming jargon as semicontinuous.  For whatever reason, the term is not well publicized; I taught a graduate course on integer programming for years without coming across it.  Some solvers can model these variables directly.  For instance, the C++ and Java APIs for CPLEX contain a class named IloSemiContVar.  I'm not sure which other solvers directly support semicontinuous variables. For solvers that do not recognize semicontinuous variables, there is a generic formulation.  You introduce a binary variable $y$ that will be 1 if $x$ is nonzero and 0 if $x$ is zero, and you add the contraints $x\ge ay$ and $x\le by$.

AMPL allows you to define a semicontinuous variable without explicitly adding constraints. Our example could be written in AMPL as -->
param a > 0;
param b > a;
var x in {0} union interval[a, b];
Currently, AMPL converts that into the pair of constraints mentioned above, even if the chosen solver is CPLEX.  Whether that will change in the future I do not know.

One potential virtue of directly modelling semicontinuous variables is that it may lead to a more parsimonious formulation.  If the solver understands semicontinuous variables, it can modify the branching logic to separate nodes based on $x=0$ versus $x\ge a$ without having to add the binary variable $y$ and the extra constraints.

Sunday, March 13, 2011

Random Sampling in AMPL

This came up on the AMPL user group; I figured I'd "double dip" and post an edited version of my solution here.  The question was how to randomly sample without replacement in AMPL.  The code below shows one way.

#
# randomly sample items from a set, with or without replacement
#
set UNIVERSE ordered;  # set of items to sample -- must be ordered
param ssize integer > 0;  # sample size desired
param sample {1..ssize} symbolic;  # will hold generated sample

data; # data below is just for demonstration purposes
set UNIVERSE := 'cat', 'dog', 'bird', 'fish', 'turtle', 'snake', 'frog', 'rock';
param ssize := 4;
model;  # back to the code

let ssize := min(ssize, card(UNIVERSE));  # omit if sampling with replacement
param item symbolic;  # dummy parameter to hold each selection
for {i in 1..ssize} {
  let item := member(ceil(Uniform(0, card(UNIVERSE))), UNIVERSE);
  let sample[i] := item;
  let UNIVERSE := UNIVERSE diff {item};  # omit if sampling with replacement
}

display sample; 

I might mention a few things about the code:
  • I used symbolic data in the example. If you are sampling from a set of numbers, you can omit the two instances of the keyword symbolic.
  • If you want to sample with replacement, there are two lines (marked by comments) that you need to omit.
  • The code puts the sample in a (vector) parameter.  If you are sampling without replacement, an alternative is to put the code into a set.  Define the set with a default value of {}, and use the union operator with second argument {item} to add the item to the sample set inside the loop.

Saturday, February 26, 2011

Finding Multiple Solutions in a Binary MIP

I'm going to post some AMPL code here that illustrates how to find more than one solution to a mixed integer program in which all integer variables are binary.  The problem I'll use is a multiple knapsack where the supply of each item is one.  The AMPL code should be pretty self-explanatory, with the possible exception of the exclusion constraints.  The notion for those is as follows.  Suppose that my binary variables are $x\in \{0,1\}^d$.  I solve the original model and obtain a solution with $x=\tilde{x}$.  I now add the following constraint:$$\sum_{i : \tilde{x}_i = 0} x_i + \sum_{i : \tilde{x}_i =1}(1-x_i)\ge 1.$$Note that this constraint forces at least one of the binary variables to change its value. (Exclusion constraints for general integer variables are harder to come by.  The only way I know off-hand involves binary expansions of the general integer variables.)

The AMPL code solves the model in a loop, recording each solution and adding an exclusion constraint.  The small example I created happens by chance to have multiple optimal solutions, enough so that the five solutions generated by the code are all optimal (but distinct).  More generally, you could get a sequence of increasingly suboptimal solutions.

Three things are perhaps noteworthy in the code.  The first is the check for infeasibility:  if you've already enumerated all possible solutions (and excluded them), the model will become infeasible.  The second is that I gave default values for the exclusion constraint coefficients (and right hand sides).  AMPL otherwise solved the base model (with no exclusion constraints) but then nagged me when I tried to record the result, because the first exclusion constraint used undefined coefficients.  Apparently recording the result (or even trying to display it) triggered a reevaluation of the model.  I think I could have ducked that by waiting to update nfound, but default values were an easy fix.  Third, note that I did not check whether a binary variable equaled 0 or 1; that might bump into rounding issues.  Saying that anything above (below) 0.5 is a 1 (0) is a very safe test.  (If your rounding error is that bad, you need to clean the beads on your abacus.)

One last comment is that this is not the most efficient way to find multiple solutions, in that it requires the (expanded) model to be solved ab initio each time.  If the solver supports callbacks, a more efficient approach is to use callbacks to record incumbents, reject them, and add exclusion constraints on the fly.  CPLEX has a "solution pool" feature that lets you track multiple solutions (and if CPLEX has it, my guess is that other commercial solvers do or soon will).  Using that is probably even  more efficient than using callbacks.  If the model is easy enough to solve, though, or if your solver lacks those features, the method demonstrated below may be a good choice.

Here is the code:
# Example of how to generate and record multiple solutions to a binary MIP.
#
# The sample problem is a 0-1 multiple knapsack.
#
set ITEMS;  # items available for the knapsacks
set KNAPSACKS;  # available knapsacks
param value {ITEMS} > 0;  # item values
param size {ITEMS} > 0;  # item sizes
param capacity {KNAPSACKS} > 0;  # knapsack capacities
var Pack {ITEMS, KNAPSACKS} binary;
  # Pack[i,k] = 1 iff item i is packed in knapsack k
maximize LoadValue: sum {i in ITEMS, k in KNAPSACKS} value[i]*Pack[i,k];
  # maximize value of load
s.t. Capacity {k in KNAPSACKS}:
  sum {i in ITEMS} size[i]*Pack[i,k]  <= capacity[k];  # capacity limits
s.t. Supply {i in ITEMS}:
  sum {k in KNAPSACKS} Pack[i,k]  <= 1;  # only one of each item available

option solver cplexamp;  # to use CPLEX as the solver

data;  # ordinarily I'd use a separate file for the data
set ITEMS := 1 2 3 4 5 6 7 8 9;
set KNAPSACKS := 1 2;
param capacity := 1 6 2 4;
param:  value  size  :=
   1       6     3
   2       5     2
   3       4     4
   4       3     1
   5       2     1
   6       4     5
   7       7     6
   8       2     3
   9       1     1  ;

model;

# we now create some structure to record multiple solutions
param want := 5;  # number of solutions we want
set SOLUTIONS := 1..want;  # index set for solutions
param packed {ITEMS, KNAPSACKS, SOLUTIONS};
  # packed[i,k,s] will be 1 if solution s packed item i in knapsack k
param payoff {SOLUTIONS};  # objective values of recorded solutions
param coef {ITEMS, KNAPSACKS, SOLUTIONS} default 0;
  # coefficients of solution exclusion constraints
param rhs {SOLUTIONS} default 0;  # right hand sides of exclusion constraints
param nfound default 0;  # number of solutions found

# we add constraints to the previous model to exclude
# solutions we've already seen
s.t. Exclude {s in 1..nfound}:
  sum {i in ITEMS, k in KNAPSACKS} coef[i,k,s]*Pack[i,k] >= rhs[s];

# solve until either the problem becomes infeasible or
# the right number of solutions has been found
for {s in SOLUTIONS} {
  solve;  # solve the model
  if solve_result = 'infeasible' then break;
    # if the model has become infeasible, give up
  let nfound := s;  # bump the counter for solutions found
  let payoff[s] := LoadValue.val;
  let rhs[s] := 1;
  # create an exclusion constraint that forces at least
  # one binary variable to change value
  for {i in ITEMS, k in KNAPSACKS} {
    if Pack[i,k] > 0.5 then {
      # treat this as a 1
      let packed[i,k,s] := 1;
      let coef[i,k,s] := -1;
      let rhs[s] := rhs[s] - 1;
    }
    else {
      # treat this as a 0
      let packed[i,k,s] := 0;
      let coef[i,k,s] := 1;
    }
  }
}