I have a problem of optimisation nearly "linear programming" where the
constraints are linear but the optimisation function is not linear :

x1 * w1 + y2 * h2 <= L
x4 * w1 + y3 * h2 <= L
x4 <= x1
y2 <= y3
y1 * h1 + y4 * h1 <= l
x2 * w2 + x3 * w2 <= l
y1 * h1 + x3 * w2 <= l
x4 * w1 + y2 * h2 <= l
All constants h1,h2,w1,w3,L,l  positive and integer
All variable positive of null and integer
Maximise x1*y1+x2*y2+x3*y3+x4*y4


A C/C++ routine or a description of an algorithm that works, ready to
be sent to a programmer, is the answer searched.

---

Clarification of Question by ocuderc-ga on 29 Jul 2006 00:59 PDT

Thank you, mapplekiwi.  Constants are from 2 to 1000, and I use now a
big loop, and it solves in a few seconds. But what I need is about
1/10 second.
I am ready to explain the origin of the problem but only privately
(send me an email at pcouderc at tol.fr)
Like the very near simplex algorithm, the set of solution is inside a
polyhedron (polytope...) defined in the space by the set of
inequations (See more explanation about in
http://en.wikipedia.org/wiki/Simplex_algorithm)
I suppose (without being sure) that we can admit that the best
solution (or  at least a good solution) can be found on one vertex of
this polyhedron.
In this case, as the number of inequations is small, the number of
vertices is fairly small too. So if I find a way to enemerate all
these vertices, I can easily (that is quickly)  compute the cost on
each of these points, and deduce the best cost.
Moreover, the set of the vertices of the polyedron is defined by a
subset of the points of intersections of the hyperplans defined by the
inequations. So the main problem is to find these intersections. As
the number of inequations is fairly small, the number of these points
will be fairly small too.

---

Clarification of Question by ocuderc-ga on 31 Jul 2006 03:58 PDT

About data, one can say that usually "w1" is equal or not too far from
"w2" or "h2", and b, respectiveley is equal or not too far from "h2"
or "w2"...
But what is "not too far" is "not too well" defined...

In your statement I assume that w3 should be replaced by w2.
Also, "of null" means "or zero".

I rewrote the problem in a more standard form:

Maximize 

  x1*y1 + x2*y2 + x3*y3 + x4*y4

subject to the constraints

  a*x1 + d*y2 <= p
  a*x4 + d*y3 <= p
  x4 <= x1
  y2 <= y3
  c*y1 + c*y4 <= q
  b*x2 + b*x3 <= q
  c*y1 + b*x3 <= q
  a*x4 + d*y2 <= q

Your problem is technically called a Quadratic Integer Progamming problem.
The web site http://www.ici.ro/camo/hint.htm refers to some packages.

How large are the constants? If they are small (say, less than 100)
then an 8-nested iteration will solve in a few seconds (most of the
time):

 for x3 = 0 .. q / b
  for y1 = 0 .. (q - b * x3) / c
    for x2 = 0 .. (q - b * x3) / b
      for y4 = 0 .. (q - c * y1) / c
        for x1 = 0 .. p / a
          for y2 = 0 .. (p - a * x1) / d
            for x4 = 0 .. min((q - d * y2) / a, x1)
              for y3 = y2 .. (p - a * x4) / d
                z = x1*y1 + x2*y2 + x3*y3 + x4*y4
                if z > zmax
                  zmax = z

I experimented with many random problems and found that many solutions
contain lots of zeros. I found only a few cases with all non-zeros,
example:

a=11 b=31 c=19 d=7 p=74 q=122 ==> x1=6 x2=1 x3=2 x4=1 y1=3 y2=1 y3=9 y4=3   

Out of curiosity can you please describe the origin of this problem.

Because x2 and y4 occur only once the 8-nested iteration is equivalent
to this 6-nested iteration:

for x1 = 0 .. p / a
  for y2 = 0 .. (p - a * x1) / d
    for x4 = 0 .. min (x1, (q - d * y2) / a)
      for y3 = y2 .. (p - a * x4) / d
        for y1 = 0 .. q / c
          for x3 = 0 .. (q - c * y1) / b                               
            x2 = q / b - x3
            y4 = q / c - y1
            z = x1 * y1 + x2 * y2 + x3 * y3 + x4 * y4
            if z > zmax
              zmax = z

Note that "/" means integer division.

A few notes:
1.  c*y1 + c*y4 <= q is equivalent to y1 + y4 <= q/c
2.  b*x2 + b*x3 <= q is equivalent to x2 + x3 <= q/b
3.  Is it true that p >= q, else a*x4 + d*y2 <= q is redundant.

Hi ocuderc

Thanks for the information.

You say that your big loop takes a few seconds.
How does this compare to my 8-nested loop? To my 6-nested loop?

Note that rounding the optimal solution of the general (non-integer)
problem is not always the optimal integer solution. See:

   www.sce.carleton.ca/faculty/chinneck/po/Chapter12.pdf

Also, unlike linear programming, the optimum of a quadratic function
is not always at a vertex.

I don't know of any efficient algorithms for listing all the integer
points close to and within the polytope. Perhaps some industrial
strength Quadratic Integer Programming systems contains such methods,
and can solve your problem using Branch & Bound.

Your problem can be solved more efficiently than the 8- or 6-nested
loops by noting that the problem can be separated into 2 parts with
disjoint variables:

Constraint 1

    a*x1 + d*y2 <= p
    a*x4 + d*y3 <= p
    a*x4 + d*y2 <= q
    x4 <= x1
    y2 <= y3

Constraint 2

    c*y1 + b*x3 <= q
    y1 + y4 <= q/c
    x2 + x3 <= q/b

Note that the objective function = x1*y1 + x4*y4+ y2*x2 + y3*x3 is
monotone increasing, and consists of independent sets of these
variables. Therefore, to find an optimal solution, first list the
integer boundaries of Constraint 1 and for each list the integer
boundaries of Constraint 2. Because Constraint 2 has single variables:
x2 and y4, you need not test !Con2(....).

Here is an improved algorithm:

for x1 = 0 ... p / a
  for y2 = 0 ... min (q / d, (p - a * x1) / d)
    for x4 = 0 ... min (x1, (p - d * y2) / a, (q - d * y2) / a)
      for y3 = y2 ... (p - a * x4) / d
        if (  !con1 (x1+1, x4,   y2,   y3    ) &
              !con1 (x1,   x4+1, y2,   y3    ) &
              !con1 (x1,   x4,   y2+1, y3    ) &
              !con1 (x1,   x4,   y2,   y3 + 1) )
           for y1 = 0 ... q / c
             for x3 = 0 ... (q - c * y1) / b
                x2 = q / b - x3
                y4 = q / c - y1
                z = x1*y1 + x2*y2 + x3*y3 + x4*y4
                if (z > zmax)
                  zmax = z
                  
This method is much faster than my previous algorithms, but sometimes
(rarely) takes longer. Example: In a=2 b=67 c=79 d=12 p=849 q=801
Constraint 1 contains only 2278 points, but takes a long time to find
them. Constraint 2 only contains 72 points.
               
Regards
maplekiwi

- As the function is monotone, I think it would be easy to prove that
the maximum is on a vertex of the polytope. Or more exactly, I think
it would be easy to prove that if the the maximum is not on a vertex,
the function is not monotone. Your use of the "monotone" word gives me
more confidence on this assertion.
- I cannot use a commercial product, because I intend to redistribute
the solution, and I have not the market to pay for it with each
licence of my software.
- I think that "Branch and Bound" is a good candidate for my problem
- I see your last proposal of algorithm with great interest. I have
not soon tried it because, the problem is very complex to translate (I
have formalised it this way, but it is more complicated) , and as the
full loop will be too big, I have anyway to find an heuristic to
eliminate some cases. So I intend to explore the ways before starting
programming.
- I would prefer to find a more generic solution, as I have
alternative problems with some variations in the constraints.

1. The following improvement yields a significant speed up.
It replaces the first 4-loop with a 3-loop, because it is necessary to
only consider the max value of y3. There may still be unnecessary
points: (x1,y2,x4,y3) so you should continue the constraint1 tests.

for x1 = 0 ... p / a
  for y2 = 0 ... min (q / d, (p - a * x1) / d)
    for x4 = 0 ... min (x1, (p - d * y2) / a, (q - d * y2) / a)
       y3 = (p - a * x4) / d; // max value of y3 <= (p-a*x4)/d
       if  !con1 (x1 + 1, x4, y2, y3) &&
           !con1 (x1, x4 + 1, y2, y3) &&
           !con1 (x1, x4, y2 + 1, y3)
           for y1 = 0 ... q / c
             for x3 = 0 ... (q - c * y1) / b
                x2 = q / b - x3
                y4 = q / c - y1
                z = x1 * y1 + x2 * y2 + x3 * y3 + x4 * y4
                ...etc...         


2. The maximum of a monotone function is not always at a vertex. Consider:

   maximize: xy
   constraints: x >= 0.  y >= 0,  x + 2y <= 4, 2x + y <= 6
   The max occurs at (2,1)
   Try this at:  pirate.shu.edu/~wachsmut/

"monotone increasing" or strictly, non-decreasing means if any
variable increases then the objective function does not decrease.

3. I appreciate your requirement for a generic solution.
   However, keep in mind that integer linear programming is NP-complete

   see:  en.wikipedia.org/wiki/Integer_programming

Currently, I use a 3 level loop by el

Currently I use a 3-level loop which seems to me very near of what you
you have designed (except that you found it mathematically, and for me
it is by imaginating the way I should do it by "hand"). Anyway, it is
still too long and I am obliged to eliminate some tests, and I use
(very) basic heuristic for that.
Anyway, I shall examine your solution, which explores all the cases,
to see if in fact I do not explore cases that cannot give results.

I do not understand well the numerical example that you give but it is
clear that for n=x*y
and a constraint like x+y<=12 and for example x<=10 and y<=10, the
maximum n will be when x=y=6 very far from the vertices....
So my assumption the the maximum of my function x1*y1+x2*y2+x3*y2+x4*y4 is false.
But maybe in this case I could explore evolutionary method such as
"ant colony".... but first integer programming should be better.

It would be interesting to compare your 3-level "heuristic" loop with
my latest "complete" algorithm. I expect that your loop may get close
to a maximum, whereas mine guarantees to find the maximum. On random
problems with constants < 1000 my algorithm produces a solution very
quickly.

You provided a simpler counterexample that the max is not always at a vertex:
maximize xy subject to x >= 0 y >= 0 x+y <= 12. (You don't need x <=
10 or y <= 10.) My example has a quadrilateral feasible region and the
max occurs in the middle of one of the slanted sides.

Given the inherent complexity (NP-complete) of integer programming
problems, heuristic methods to solve general quadratic integer
programming problems may have to suffice to obtain a good (but not
optimal) solution in a reasonable time.

My current loop with my basic(!) heuristic is something like :
for x1 = 0 ... min(p / a, MAXALLOWED1)
  for y2 = 0 ... min (q / d, (p - a * x1) / d, MAXALLOWED2)
    for x4 = 0 ... min (x1, (p - d * y2) / a, (q - d * y2) / a, MAXALLOWED3)
       select 1 value for other variables...

Thank you for brain storming me so, now I go to see "branch and bound"
in wikipedia...