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
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:
- i is an item with color blue
- j is an item with color not blue
- 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
- How to construct a complex constraint in PuLP Python, https://stackoverflow.com/questions/58770011/how-to-construct-a-complex-constraint-in-pulp-python
Or just solve three separate problems (in parallel) where the blue shape is fixed.
ReplyDelete