Monday, November 11, 2019

Interesting constraint

Problem statement


In [1] the following problem is proposed:

I'm trying to solve a knapsack-style optimization problem with additional complexity.
Here is a simple example. I'm trying to select 5 items that maximize value. 2 of the items must be orange, one must be blue, one must be yellow, and one must be red. This is straightforward. However, I want to add a constraint that the selected yellow, red, and orange items can only have one shape in common with the selected blue item.

Data


The example data looks like:


 item   color     shape  value
    A    blue    circle   0.454
    B  yellow    square   0.570
    C     red  triangle   0.789
    D     red    circle   0.718
    E     red    square   0.828
    F  orange    square   0.709
    G    blue    circle   0.696
    H  orange    square   0.285
    I  orange    square   0.698
    J  orange  triangle   0.861
    K    blue  triangle   0.658
    L  yellow    circle   0.819
    M    blue    square   0.352
    N  orange    circle   0.883
    O  yellow  triangle   0.755


The complicated constraint is a bit difficult to grasp right away. The two solutions listed below may help to understand how I interpreted this.

Derived data


Let's see if we can model this. First we slice and dice the data a bit to make the modeling a bit easier. Here is some derived data:


----     65 SET i  item

A,    B,    C,    D,    E,    F,    G,    H,    I,    J,    K,    L,    M,    N,    O


----     65 SET c  color

blue  ,    yellow,    red   ,    orange


----     65 SET s  shape

circle  ,    square  ,    triangle


----     65 SET ICS  (i,c,s) combinations

              circle      square    triangle

A.blue           YES
B.yellow                     YES
C.red                                    YES
D.red            YES
E.red                        YES
F.orange                     YES
G.blue           YES
H.orange                     YES
I.orange                     YES
J.orange                                 YES
K.blue                                   YES
L.yellow         YES
M.blue                       YES
N.orange         YES
O.yellow                                 YES


----     65 SET IC  (i,c) combinations

         blue      yellow         red      orange

A         YES
B                     YES
C                                 YES
D                                 YES
E                                 YES
F                                             YES
G         YES
H                                             YES
I                                             YES
J                                             YES
K         YES
L                     YES
M         YES
N                                             YES
O                     YES


----     65 SET IS  (i,s) combinations

       circle      square    triangle

A         YES
B                     YES
C                                 YES
D         YES
E                     YES
F                     YES
G         YES
H                     YES
I                     YES
J                                 YES
K                                 YES
L         YES
M                     YES
N         YES
O                                 YES


----     65 SET blue  blue colored items

A,    G,    K,    M


----     65 PARAMETER value  value of item

A 0.454,    B 0.570,    C 0.789,    D 0.718,    E 0.828,    F 0.709,    G 0.696,    H 0.285,    I 0.698,    J 0.861
K 0.658,    L 0.819,    M 0.352,    N 0.883,    O 0.755


I often introduce data in more convenient formats. This will help us in writing the model in a more readable way. We can develop and debug this derived data before we work on the model equations itself.

Note that the set CS(c,s) is complete. However, I will assume that there is a possibility that this set has some missing entries. In other words, I will not assume that all combinations of colors and shapes exist in the data.

High-level model I


Let's introduce the following zero-one variables:\[\begin{align} & x_i = \begin{cases} 1 & \text{if item $i$ is selected}\\ 0 & \text{otherwise}\end{cases} \\ & y_{c,s} = \begin{cases} 1 & \text{if items with color $c$ and shape $s$ are selected}\\ 0 & \text{otherwise} \end{cases}\end{align}\]


My high-level model is:

High-level Model 1
\[\begin{align}\max & \sum_i \color{darkblue}{\mathit{Value}}_i \cdot \color{darkred}x_i \\ &\sum_i \color{darkred}x_i = \color{darkblue}{\mathit{NumItems}}\\ &\sum_{i | \color{darkblue}{\mathit{IC}}(i,c)} \color{darkred}x_i = \color{darkblue}{\mathit{NumColor}}_c && \forall c\\ & \color{darkred}y_{c,s} = \max_{i|\color{darkblue}{\mathit{ICS}}(i,c,s)} \color{darkred}x_i && \forall c,s|\color{darkblue}{\mathit{CS}}(c,s)\\  & \color{darkred}y_{\color{darkblue}{\mathit{blue}},s} = 1 \Rightarrow \sum_{c|\color{darkblue}{\mathit{YRO}}(c)} \color{darkred}y_{c,s} \le 1 && \forall s \\&\color{darkred}x_i, \color{darkred}y_{c,s} \in \{0,1\} \end{align}\]

The max constraint implements the definition of the \(y_{c,s}\) variables: if any of the selected items has color/shape combination \((c,s)\)  then \(y_{c,s}=1\) (and else it stays zero). The implication constraint says: if we have a blue shape \(s\), then there can be only one shape of type \(s\) of other color. The model is a bit complicated with all subsets and indexing because I wanted to be precise. No hand-waving. This helps when implementing it.

MIP Model I


This is not yet a MIP model, but translation of the above model into a normal MIP is not too difficult.


Mixed Integer Programming Model 1
\[\begin{align}\max & \sum_i \color{darkblue}{\mathit{Value}}_i \cdot \color{darkred}x_i \\ &\sum_i \color{darkred}x_i = \color{darkblue}{\mathit{NumItems}}\\ &\sum_{i | \color{darkblue}{\mathit{IC}}(i,c)} \color{darkred}x_i = \color{darkblue}{\mathit{NumColor}}_c && \forall c\\ & \color{darkred}y_{c,s} \ge \color{darkred}x_i && \forall i,c,s|\color{darkblue}{\mathit{ICS}}(i,c,s)\\  & \sum_{c|\color{darkblue}{\mathit{YRO}}(c)} \color{darkred}y_{c,s} \le 1 + \color{darkblue}M(1-\color{darkred}y_{\color{darkblue}{\mathit{blue}},s}) && \forall s \\&\color{darkred}x_i, \color{darkred}y_{c,s} \in \{0,1\} \end{align}\]

The max constraint has been replaced by a single inequality. This works in this case as we are only interested in \(y\)'s that are forced to be one. Those guys will limit blue shape constraint. This type of bounding is quite common, but it always requires a bit of thinking through in order to establish if this is allowed for the specific case. The blue shape implication constraint itself is rewritten as a big-M inequality. A good value for \(M\) is not very difficult to establish: the number of colors minus blue minus 1.


MIP Model 2


We can actually bypass the variable \(y\) and directly formulate the limit on selected shapes in terms of variables \(x\). For this we develop a set SameShape(i,j) with:

  1. i is an item with color blue
  2. j is an item with color not blue
  3. i and j have the same shape

The set SameShape looks like:


----     71 SET SameShape  blue i and non-blue j with same shape

            B           C           D           E           F           H           I           J           L

A                                 YES                                                                     YES
G                                 YES                                                                     YES
K                     YES                                                                     YES
M         YES                                 YES         YES         YES         YES

+           N           O

A         YES
G         YES
K                     YES


This means for the first row: if item A (a blue circle) is considered, then other non-blue circles are items D, L and N. With this we can write:


Mixed Integer Programming Model 2
\[\begin{align}\max & \sum_i \color{darkblue}{\mathit{Value}}_i \cdot \color{darkred}x_i \\ &\sum_i \color{darkred}x_i = \color{darkblue}{\mathit{NumItems}}\\ &\sum_{i | \color{darkblue}{\mathit{IC}}(i,c)} \color{darkred}x_i = \color{darkblue}{\mathit{NumColor}}_c && \forall c\\   &  \sum_{j|\color{darkblue}{\mathit{SameShape}}(i,j)} \color{darkred}x_j \le 1 + \color{darkblue}M(1-\color{darkred}x_i) && \forall i | \color{darkblue}{\mathit{Blue}}_i  \\&\color{darkred}x_i \in \{0,1\} \end{align}\]

Solution without "blue shape" constraint


For comparison, let's first run the model without the complex "blue" shape constraints (the last two constraints). That gives:


----     90 VARIABLE z.L                   =        4.087  obj

----     90 VARIABLE x.L  select item

E 1.000,    G 1.000,    J 1.000,    L 1.000,    N 1.000


----     90 SET selected  

              circle      square    triangle

E.red                        YES
G.blue           YES
J.orange                                 YES
L.yellow         YES
N.orange         YES

We see that the blue shape (circle) is also selected as orange and yellow.

Solution with "blue shape" constraint


With the additional constraints, we see:


----     95 VARIABLE z.L                   =        4.049  obj

----     95 VARIABLE x.L  select item

E 1.000,    J 1.000,    K 1.000,    L 1.000,    N 1.000


----     95 SET selected  

              circle      square    triangle

E.red                        YES
J.orange                                 YES
K.blue                                   YES
L.yellow         YES
N.orange         YES

Now we see the blue shape is a triangle. We only have another triangle of color orange.

Obviously, as a sanity check, the objective reduces a bit when we introduce this constraint.


Appendix: Model formulated in GAMS


The complete GAMS model looks like:

$ontext

  
I'm trying to solve a napsack-style optimization problem with additional complexity.

  
Here is a simple example. I'm trying to select 5 items that maximize value. 2 of the
  
items must be orange, one must be blue, one must be yellow, and one must be red.
  
This is straightforward. However, I want to add a constraint that the selected yellow,
  
red, and orange items can only have one shape in common with the selected blue item.

$offtext


set
   i
'item' /A*O/
   c
'color' /blue,yellow,red,orange/
   s
'shape' /circle,square,triangle/
;

parameters
    data(i,c,s)
'value' /
        
A . blue   .  circle     0.454
        
B . yellow .  square     0.570
        
C . red    .  triangle   0.789
        
D . red    .  circle     0.718
        
E . red    .  square     0.828
        
F . orange .  square     0.709
        
G . blue   .  circle     0.696
        
H . orange .  square     0.285
        
I . orange .  square     0.698
        
J . orange .  triangle   0.861
        
K . blue   .  triangle   0.658
        
L . yellow .  circle     0.819
        
M . blue   .  square     0.352
        
N . orange .  circle     0.883
        
O . yellow .  triangle   0.755
        
/
    NumItems
'number of items to select' /5/
    NumColor(c)
'required number of each color' /
        
orange 2
        
red    1
        
blue   1
        
yellow 1
       
/
;

alias(i,j);

sets
  ICS(i,c,s)
'(i,c,s) combinations'
  IC(i,c)   
'(i,c) combinations'
  IS(i,s)   
'(i,s) combinations'
  blue(i)   
'blue colored items'
;
parameter
  value(i)  
'value of item'
  M         
'big-M'
;

ICS(i,c,s) = data(i,c,s);
IC(i,c) =
sum(ICS(i,c,s),1);
IS(i,s) =
sum(ICS(i,c,s),1);
blue(i) = IC(i,
"blue");
value(i) =
sum((c,s),data(i,c,s));
display i,c,s,ICS,IC,IS,blue,value;

M =
card(s)-2;

set SameShape(i,j) 'blue i and non-blue j with same shape';
SameShape(i,j)$(blue(i)
and not blue(j)) = sum(s, IS(i,s) and IS(j,s));
display SameShape;

binary variable x(i) 'select item';
variable z 'obj';

equations
   obj      
'objective'
   count    
'count number of selected items'
   countcolor(c) 
'count selected items for each color'
   limit(i)
;
obj.. z =e=
sum(i, value(i)*x(i));
count..
sum(i, x(i)) =e= numitems;
countcolor(c)..
sum(IC(i,c), x(i)) =e= numcolor(c);
limit(i)$blue(i)..
sum(sameshape(i,j), x(j)) =l= 1 + M*(1-x(i));

option optcr=0;
set selected(i,c,s) 'reporting';

model m1 /obj,count,countcolor/;
solve m1 maximizing z using mip;
selected(i,c,s)$ICS(i,c,s) = x.l(i)>0.5;
display z.l,x.l,selected;

model m2 /all/;
solve m2 maximizing z using mip;
selected(i,c,s)$ICS(i,c,s) = x.l(i)>0.5;
display z.l,x.l,selected;



A second exercise would be to write this in Python using PuLP.


Appendix: PuLP model


import pulp
import io
import pandas

# -------- data ----------------

df = pandas.read_csv(io.StringIO(
'''
item,color,shape,value
A,blue,circle,0.454
B,yellow,square,0.570
C,red,triangle,0.789
D,red,circle,0.718
E,red,square,0.828
F,orange,square,0.709
G,blue,circle,0.696
H,orange,square,0.285
I,orange,square,0.698
J,orange,triangle,0.861
K,blue,triangle,0.658
L,yellow,circle,0.819
M,blue,square,0.352
N,orange,circle,0.883
O,yellow,triangle,0.755
'''))
print(df)

numItems =
5
numColor = {
'orange':2, 'red':1, 'blue':1, 'yellow':1}

# --------- derived data ------------

Items = df[
"item"]
Colors = df[
"color"].unique()
Shapes = df[
"shape"].unique()
YRO = Colors[Colors !=
'blue']

value = dict(zip(Items,df[
"value"]))
color = dict(zip(Items,df[
"color"]))
shape = dict(zip(Items,df[
"shape"]))

M = Shapes.size -
2

SameShape = {
   (i,j)
for i in Items for j in Items
        
if (color[i]=='blue') and (color[j]!='blue') and
             (shape[j]==shape[i]) }
print('SameShape:',SameShape)

# ----------- Model -----------

p = pulp.LpProblem(
"MIP problem", pulp.LpMaximize)
x = pulp.LpVariable.dicts(
"x",Items,0,1,pulp.LpInteger)

p += pulp.lpSum([x[i]*value[i]
for i in Items])
p += pulp.lpSum([x[i]
for i in Items]) == numItems

for c in Colors:
    p += pulp.lpSum([x[i]
for i in Items if color[i]==c]) == numColor[c]

for i in Items:
   
if color[i]=="blue":
         p += pulp.lpSum([x[j]
for j in Items if (i,j) in SameShape]) <= 1 + M*(1-x[i])

p.solve(pulp.PULP_CBC_CMD(msg=
1))

#----------- results --------------

print("Status:", pulp.LpStatus[p.status])
print("Objective:",pulp.value(p.objective))

selected = []
for i in Items:
   
if pulp.value(x[i]) > 0.5:
         selected += [[i, color[i], shape[i]]]
selected = pandas.DataFrame(selected, columns=[
'item','color','shape'])
print(selected)



Note that SameShape is a Python set. We could have dropped SameShape and implemented the logic of this set directly in the blue shape constraint. However I believe it is often a good thing to keep constraints as simple as possible. They are more difficult to debug. The set SameShape can easily be inspected by a print statement, even before we have written the rest of the MIP model.

The output looks like:


   item   color     shape  value
0     A    blue    circle  0.454
1     B  yellow    square  0.570
2     C     red  triangle  0.789
3     D     red    circle  0.718
4     E     red    square  0.828
5     F  orange    square  0.709
6     G    blue    circle  0.696
7     H  orange    square  0.285
8     I  orange    square  0.698
9     J  orange  triangle  0.861
10    K    blue  triangle  0.658
11    L  yellow    circle  0.819
12    M    blue    square  0.352
13    N  orange    circle  0.883
14    O  yellow  triangle  0.755
Welcome to the CBC MILP Solver
Version: 2.9.0
Build Date: Feb 12 2015

command line - D:\Python\Python37\lib\site-packages\pulp\solverdir\cbc\win\64\cbc.exe 3c0bdd177daf444ba924d26442c171ba-pulp.mps max branch printingOptions all solution 3c0bdd177daf444ba924d26442c171ba-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 14 COLUMNS
At line 108 RHS
At line 118 BOUNDS
At line 134 ENDATA
Problem MODEL has 9 rows, 15 columns and 48 elements
Coin0008I MODEL read with 0 errors
Continuous objective value is 4.049 - 0.00 seconds
Cgl0004I processed model has 9 rows, 15 columns (15 integer (15 of which binary)) and 48 elements
Cutoff increment increased from 1e-005 to 0.000999
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of -4.049
Cbc0038I Before mini branch and bound, 15 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.01 seconds)
Cbc0038I After 0.01 seconds - Feasibility pump exiting with objective of -4.049 - took 0.00 seconds
Cbc0012I Integer solution of -4.049 found by feasibility pump after 0 iterations and 0 nodes (0.01 seconds)
Cbc0001I Search completed - best objective -4.049, took 0 iterations and 0 nodes (0.01 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from -4.049 to -4.049
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                4.04900000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.02
Time (Wallclock seconds):       0.02

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.07   (Wallclock seconds):       0.07

Status: Optimal
Objective: 4.0489999999999995
  item   color     shape
0    E     red    square
1    J  orange  triangle
2    K    blue  triangle
3    L  yellow    circle
4    N  orange    circle



Some notes:

  • This does not look so bad. The PuLP model follows the mathematical model quite closely.
  • The main differences with the GAMS model: we need for loops here and there (GAMS has implicit loops), and the variable y is dense. I don't think PuLP has direct facilities for sparse variables. In this case this is not a problem.
  • The CBC version that comes with PuLP is rather old.
  • Having said that, CBC is doing pretty good on this model: 0 nodes and 0 simplex iterations.

References


1 comment:

  1. Or just solve three separate problems (in parallel) where the blue shape is fixed.

    ReplyDelete