Tuesday, March 31, 2009

MTSP/VRP example

> You have an example of a VRP (Vehicle Routing Problem)?

This is an example of a Multiple traveling salesman problem or uncapacitated VRP. A very simple formulation though.

 

$ontext

  
bayg29 29 city problem  with five salesman

  
Erwin Kalvelagen, Amsterdam Optimization

  
References:
      
Bektas, T. (2006). The multiple traveling salesman problem: an overview of formulations and
       
solution procedures. OMEGA: The International Journal of Management Science, 34(3), 209-219.

$offtext


$set startcity i13


scalar
   m
'number of salesman' /5/
   p
'number of cities to visit' /6/
;

sets
   i
/i1*i29/
;
alias (i,j);

table c(i,j)
'distance matrix (KM)'


     
i1  i2  i3  i4  i5  i6  i7  i8  i9 i10 i11 i12 i13 i14 i15 i16 i17 i18 i19 i20 i21 i22 i23 i24 i25 i26 i27 i28 i29
 
i1       97 205 139  86  60 220  65 111 115 227  95  82 225 168 103 266 205 149 120  58 257 152  52 180 136  82  34 145
 
i2          129 103  71 105 258 154 112  65 204 150  87 176 137 142 204 148 148  49  41 211 226 116 197  89 153 124  74
 
i3              219 125 175 386 269 134 184 313 201 215 267 248 271 274 236 272 160 151 300 350 239 322  78 276 220  60
 
i4                  167 182 180 162 208  39 102 227  60  86  34  96 129  69  58  60 120 119 192 114 110 192 136 173 173
 
i5                       51 296 150  42 131 268  88 131 245 201 175 275 218 202 119  50 281 238 131 244  51 166  95  69
 
i6                          279 114  56 150 278  46 133 266 214 162 302 242 203 146  67 300 205 111 238  98 139  52 120
 
i7                              178 328 206 147 308 172 203 165 121 251 216 122 231 249 209 111 169  72 338 144 237 331
 
i8                                  169 151 227 133 104 242 182  84 290 230 146 165 121 270  91  48 158 200  39  64 210
 
i9                                      172 309  68 169 286 242 208 315 259 240 160  90 322 260 160 281  57 192 107  90
 
i10                                         140 195  51 117  72 104 153  93  88  25  85 152 200 104 139 154 134 149 135
 
i11                                             320 146  64  68 143 106  88  81 159 219  63 216 187  88 293 191 258 272
 
i12                                                 174 311 258 196 347 288 243 192 113 345 222 144 274 124 165  71 153
 
i13                                                     144  86  57 189 128  71  71  82 176 150  56 114 168  83 115 160
 
i14                                                          61 165  51  32 105 127 201  36 254 196 136 260 212 258 234
 
i15                                                             106 110  56  49  91 153  91 197 136  94 225 151 201 205
 
i16                                                                 215 159  64 126 128 190  98  53  78 218  48 127 214
 
i17                                                                      61 155 157 235  47 305 243 186 282 261 300 252
 
i18                                                                         105 100 176  66 253 183 146 231 203 239 204
 
i19                                                                             113 152 127 150 106  52 235 112 179 221
 
i20                                                                                  79 163 220 119 164 135 152 153 114
 
i21                                                                                     236 201  90 195  90 127  84  91
 
i22                                                                                         273 226 148 296 238 291 269
 
i23                                                                                             112 130 286  74 155 291
 
i24                                                                                                 130 178  38  75 180
 
i25                                                                                                     281 120 205 270
 
i26                                                                                                         213 145  36
 
i27                                                                                                              94 217
 
i28                                                                                                                 162
 
i29
;



c(i,j) = max(c(i,j),c(j,i));

set arcs(i,j);
arcs(i,j)$c(i,j) =
yes
;

set i0(i) /%startcity%/
;

set
i2(i);
i2(i)$(
not i0(i)) = yes
;

scalar n 'number of nodes'
;
n =
card
(i);

binary variables
x(i,j);
positive variables
u(i);
variable
z;

equations

   start
   end
   assign1(i)
   assign2(j)
   sec(i,j)   
'subtour elimination'
   cost
;



start..
sum(i2(j), x('%startcity%',j)) =e= m;

end..
sum(i2(i), x(i,'%startcity%'
)) =e= m;

assign1(i2(i))..
sum
(arcs(i,j), x(i,j)) =e= 1;

assign2(i2(j))..
sum
(arcs(i,j), x(i,j)) =e= 1;

sec(arcs(i,j))$(i2(i)
and
i2(j)).. u(i) - u(j) + p*x(i,j) =L= p-1;

cost.. z =e=
sum
(arcs, c(arcs)*x(arcs));

option
optcr=0;
option
iterlim=10000000;
option
reslim=3600;

* use parallel MIP

option threads=8;
option
mip=cplex;

model mtsp/all/
;
solve
mtsp minimizing z using mip;

display
x.l;

table
xy(i,*)
             
x       y

  
i1    1150.0  1760.0
  
i2     630.0  1660.0
  
i3      40.0  2090.0
  
i4     750.0  1100.0
  
i5     750.0  2030.0
  
i6    1030.0  2070.0
  
i7    1650.0   650.0
  
i8    1490.0  1630.0
  
i9     790.0  2260.0
 
i10     710.0  1310.0
 
i11     840.0   550.0
 
i12    1170.0  2300.0
 
i13     970.0  1340.0
 
i14     510.0   700.0
 
i15     750.0   900.0
 
i16    1280.0  1200.0
 
i17     230.0   590.0
 
i18     460.0   860.0
 
i19    1040.0   950.0
 
i20     590.0  1390.0
 
i21     830.0  1770.0
 
i22     490.0   500.0
 
i23    1840.0  1240.0
 
i24    1260.0  1500.0
 
i25    1280.0   790.0
 
i26     490.0  2130.0
 
i27    1460.0  1420.0
 
i28    1260.0  1910.0
 
i29     360.0  1980.0
;



parameter tour(*,*);
set kall/k1*k10000/
;
set kk(kall) /k1/
;

loop
((i,j)$(x.l(i,j)>0.5),
   tour(kk,
'x0') = xy(i,'x'
);
   tour(kk,
'y0') = xy(i,'y'
);
   tour(kk,
'x1') = xy(j,'x'
);
   tour(kk,
'y1') = xy(j,'y'
);
   kk(kall)=kk(kall-1);
);

display
tour;



 

AMPL, OML version of GAMS Scheduling Model

In the good tradition of letting blogs link to each other, here is a blog posting with an AMPL and OML version of this GAMS model. Those are nicely formatted listings.

Sunday, March 29, 2009

Efficient Portfolio Frontier in MSF (2)

Here is an interesting but obscure reformulation that sometimes performs better than the traditional x’Qx term in a mean-variance portfolio model. We notice that the covariances can be calculated as:

covar where μi is the mean return and di,t are the historical data. We now can write:

covar2wherecovar3

This is especially a good formulation if there are more instruments than time periods (in large practical problems this is often the case as older historical data is not considered as having much predictive value). The corresponding OML model, where we solve K sub-problems at the same time, can look like (see post for more info):

Model[
   Parameters[Sets,I,K,T],
   Parameters[Integers,Selected[I]],
   Parameters[Reals,Return[I]],
   Parameters[Reals,Lambda[K]],
   Parameters[Reals,MeanAdjRet[I,T]],
   Parameters[Integers,TT[]],

   Decisions[Reals[0,Infinity],Alloc[K,I]],
   Decisions[Reals,w[K,T]],

   Constraints[
       Foreach[{k,K},FilteredSum[{i,I},Selected[i]>0,Alloc[k,i]]==1],
       Foreach[{k,K},{t,T},w[k,t]==FilteredSum[{i,I},Selected[i]>0,MeanAdjRet[i,t]*Alloc[k,i]]]
   ],

   Goals[Minimize["Overall"->
         Sum[{k,K},
              Sum[{t,T},w[k,t]^2]/TT[]
              - Lambda[k]*FilteredSum[{i,I},Selected[i]>0,Return[i]*Alloc[k,i]]
           ]
   ]]
]

As an aside: this formulation also shows how this can be linearized: replace ∑w2 by ∑|w|.

Saturday, March 28, 2009

Yet another big-M problem, don’t blame Cplex

Hi,

I posted this message at the ILOG CPLEX Forum but I guess it could be
of interest here. I know that some people watch both forums, so I
apologize for the cross-post.

Any idea is appreciated. Could it be a problem on our cplex
installation ? We are using Ubuntu (one of the 2008 versions, I can
check exactly which one if anyone finds that information important).
---

Dear all,

I am really puzzled with the following:

I have an IP model. ( …) . If I solve this model (I am using Cplex 11.2, interative
optimizer, under Ubuntu ), I obtain an "optimal" integer solution with
cost 23 (…
).

HOWEVER, I know of a solution with cost 21 ( … ). Indeed, if I
pass this solution to CPLEX as an mst file, Cplex starts using it as
the "Best Integer found" ( …
).

Am I missing anything ? Any help is appreciated,

It was suggested by another poster to use some clever but perfectly sensible options:

I verified this with CPLEX 11.0 (also on Ubuntu).  There are two ways (that I found) to get CPLEX to produce a correct solution (objective = 21) without the hot start:  turn off the presolver; or crank up the simplex feasibility tolerance (1e-9 worked for me).  Turning off the presolver seems to be the faster option.  (Disclaimer:  In both instances, I aborted the run after it found an integer-feasible solution with value 21 but before it exhausted the search tree, so it's possible one or the other option might have led to problems down the road.)

Cplex support in this forum adds more suggestions for even more exotic options (again they make perfect sense). However, I would approach this a little bit differently. First this problem is not purely a Cplex problem. Other solvers show the same problem. The reason lies actually in the formulation, which causes problems for relaxation-based branch & bound methods. It is always difficult to recognize the model structure from an LP file (please, use a proper modeling language) but there are big-M constants in the model with value 100,000. After I replaced all these instances by 10,000 Cplex shows the expected behavior: it does not stop with obj=23 but goes further to obj=21 on Ubuntu. So rather than tinkering with tolerances I would work on the model formulation and provide big-M values that are as small as possible. These values are often related to some upper bound. Experienced modelers (you can sometimes recognize them by their grey haires and thick glasses) will try to prevent using just a big arbitrary value like 1e5 but rather work on deducing a tight bound. Once you have proper formulation you will see that the model will behave better for any MIP solver. I would argue that a solver independent fix is better than using a solver specific option that works around the problem. (If you absolutely cannot provide a good value for your big-M’s you may want to look into Cplex’s indicator constraints).

Efficient Portfolio Frontier in MSF

I received a few questions about how to calculate an efficient portfolio frontier in MS Solver Foundation. The simplest algorithm is to solve a series of QP’s. (A better way may be a parametric QP analysis). As we don’t have looping constructs in OML, it is not immediately obvious how to do this with the Excel plug-in. With the API’s we could do the looping in C# and the Quadratic Programming Problem itself in OML. Here is an alternative approach: solve a bigger problem that has all the individual QP’s embedded. Conceptually this would look like:

Obj: x1’Qx1-λ1 μ'x1 +x2’Qx2-λ2 μ'x2 +x3’Qx3-λ3 μ'x3 +x4’Qx4-λ4 μ'x4
S.t. ix1,i=1      
    ix2,i=1    
      ix3,i=1  
        ix4,i=1

Note that xk is a vector rather than a scalar, and that I ignored the non-negativity conditions xk≥0. A similar approach can be used for any problem where the individual problems are relatively small. A good example is DEA where we achieve better performance by batching small problems into bigger ones.

The OML formulation of this combined portfolio model can look like:

Model[
   Parameters[Sets,I,K],
   Parameters[Integers,Selected[I]],
   Parameters[Reals,Return[I]],
   Parameters[Reals,Lambda[K]],
   Parameters[Reals,Covar[I,I]],

   Decisions[Reals[0,Infinity],Alloc[K,I]],
   Constraints[
       Foreach[{k,K},FilteredSum[{i,I},Selected[i]>0,Alloc[k,i]]==1]
   ],

   Goals[Minimize["Overall"->
         Sum[{k,K},
              FilteredSum[{i,I},Selected[i]>0,
                   FilteredSum[{j,I},Selected[j]>0,Alloc[k,i]*Covar[i,j]*Alloc[k,j]]]
              - Lambda[k]*FilteredSum[{i,I},Selected[i]>0,Return[i]*Alloc[k,i]]
           ]
   ]]
]

Here the parameter Selected indicates whether an instrument can be part of the portfolio.

Thursday, March 26, 2009

Drawing TSP Tour in GAMS

 Can I draw a TSP tour in GAMS?

Yes, use 2D vectors. Here are a few examples: a 29 city problem and a 51 city problem solved as a MIP.



Wednesday, March 25, 2009

A machine scheduling problem

Hello

I'm new to AMPL.  I have  a scheduling problem that is giving me
trouble. I have a set of unrelated machines.  I have a set of jobs.
Each job has a set of operations to perform in order.  I want to
minimize the time of completion of the last job/operation pair
processed. Below is my .mod file. In the data section I have two test
cases - one is commented and one is not.

The commented test has no problems. The solver finds the solution. The
uncommented test case gives the following:

CPLEX Error  1217: No solution exists.
CPLEX Error  1217: No solution exists.
CPLEX Error  1262: No basis exists.
CPLEX 7.1.0: integer infeasible.
1894 MIP simplex iterations
39 branch-and-bound nodes; no basis.

Like I said I am new to this type of programming. Does anybody smarter
than me see a problem with any of my constraints? Or a comment on what
usually causes error 1217 and 1262?

#################################### .mod file
set job;
set oper;
set mach;
param first > 0;
param last  > first;
set T := first .. last;

var max_val;
var asgn { job, oper, mach, T } binary;
var sched  { mach, T } >= 0;
var delay  { mach, T } >= 0;

param time { job, oper, mach } >= 0;

minimize z: max_val;

# the last scheduled task for each machine is less than or equal to
our objective
subject to max_val_constraint { m in mach }:  sched[m,last] <=
max_val;

# sched[m,t] is the completion time of the task machine m is doing at
time t
# sched[m,t] = sched[m,t] + delay[m,t] + { time to process assigned
task}
subject to mach_sched_0 { m in mach, t in T: t>1}:
  sched[m,t] = sched[m,t-1] + delay[m,t] +
    sum{ i in job, j in oper } asgn[i,j,m,t] * time[i,j,m];

subject to mach_sched_1 { m in mach }:
  sched[m,1] = delay[m,1] +
    sum{ i in job, j in oper } asgn[i,j,m,1] * time[i,j,m];

# for each job/operation pair, only assign it to be processed once
subject to jobtask_assign_cnt { i in job, j in oper }:
  sum{ m in mach, t in T } asgn[i,j,m,t] = 1;

# for each machine at time t, assign only one task to process
subject to machtime_assign_cnt { m in mach, t in T }:
  sum{ i in job, j in oper } asgn[i,j,m,t] <= 1;

# for job i oper j, make sure job i, oper j-1 has been previously assigned
subject to jobtask_ordering { i in job, j in oper, t in T: j>1 }:
  sum{ m in mach } asgn[i,j,m,t] >= sum{ m in mach } asgn[i,j-1,m,t];

subject to delay_constraint { m in mach, t in T }:
  delay[m,t] >= 0;

data;

set job  := 1 2 3 ;
set oper := 1 2 3 4;
set mach := 1 2 ;

param first := 1;
param last  := 12;

param time :=
[1,*,*]: 1 2 :=
1 1 4
2 2 2
3 3 3
4 4 4

[2,*,*]: 1 2 :=
1 2 4
2 2 1
3 1 6
4 1 2

[3,*,*]: 1 2  :=
1 3 4
2 2 3
3 3 1
4 4 3;

# test that works
#set job  := 1 2;
#set oper := 1 2;
#set mach := 1 2 3;
#set T := 1 2 3 4;
#
#param time :=
#[1,*,*]: 1 2 3 :=
#1 1 2 3
#2 1 2 3
#
#[2,*,*]: 1 2 3 :=
#1 1 2 3
#2 1 2 3;

The constraint

# for job i oper j, make sure job i, oper j-1 has been previously assigned
subject to jobtask_ordering { i in job, j in oper, t in T: j>1 }:
  sum{ m in mach } asgn[i,j,m,t] >= sum{ m in mach } asgn[i,j-1,m,t];

does not seem to be formulated correctly. It says that all operations have to be executed in the same time period t. For this data set you have four operations and only two machines available. The smaller data set had two operations which can be allocated to the available two machines. If you really want (i,j-1) to be finished before (i,j) you probably need to keep track when each task (i,j) finishes. I believe your current formulation mixes up time slots t and completion times modeled via sched[m,t]. If your execution times are small integer, then map everything to time slots. If your execution times are floating point numbers, you should look into a continuous time formulation (e.g. through x[i,j,ii,jj,m] = 1 iff (i,j) precedes (ii,jj) on machine m).

MS Solver Foundation job

It was brought to my attention that already there are (short term) Microsoft Foundation jobs available: http://www.careerbuilder.com/JobSeeker/Jobs/JobDetails.aspx?job_did=J3G7ML643DJNGKK4315. That means people start to do serious work with this!

rminlp’s

Hi all,
I am working in pharma and here I need you help for a problem related to a pharmacokinetic equation.
I would like to have the (X,Y,S, T, U, V, W) possible values that will lead to a value of Z=5.7.
Z= 1/X*1/(Y/(1+(U×V)/(W×(U+S)))+(1-Y)) )
knwing that:
0.57<X<0.79
0.00016<W<0.00032
0.89<Y<0.94
0.0936<U<0.967
V =0.05 ou 0.072
S = 29 ou 5.3

This GAMS model calculates some solutions:

variables z,x,y,u,v,w,s;
equation e;

e.. z =e= (1/X)*1/(Y/(1+(U*V)/(W*(U+S)))+(1-Y));

x.lo = 0.57;    x.up = 0.79;
w.lo = 0.00016; w.up = 0.00032;
y.lo = 0.89;    y.up = 0.94;
u.lo = 0.0936;  u.up = 0.0967;

binary variable b1,b2;
equation e1,e2;
e1.. v =e= 0.05*b1 + 0.072*(1-b1);
e2.. s =e= 29*b2 + 5.3*(1-b2);

z.fx = 5.7;

model m /all/;

b1.fx=0; b2.fx=0; solve m minimizing z using rminlp;
b1.fx=0; b2.fx=1; solve m minimizing z using rminlp;
b1.fx=1; b2.fx=0; solve m minimizing z using rminlp;
b1.fx=1; b2.fx=1; solve m minimizing z using rminlp;

I could only find feasible solutions for b2=0 (i.e. s=5.3). This was confirmed by the global solver Baron. Here they are:

                           LOWER          LEVEL          UPPER         MARGINAL

---- VAR z                  5.7000         5.7000         5.7000         1.0000     
---- VAR x                  0.5700         0.6019         0.7900          .         
---- VAR y                  0.8900         0.8900         0.9400         EPS        
---- VAR u                  0.0936         0.0936         0.0967         EPS        
---- VAR v                 -INF            0.0720        +INF             .         
---- VAR w                  0.0002         0.0003         0.0003         EPS        
---- VAR s                 -INF            5.3000        +INF             .         
---- VAR b1                  .              .              .             EPS        
---- VAR b2                  .              .              .             EPS
        

                           LOWER          LEVEL          UPPER         MARGINAL

---- VAR z                  5.7000         5.7000         5.7000         1.0000     
---- VAR x                  0.5700         0.5700         0.7900         EPS        
---- VAR y                  0.8900         0.9400         0.9400         EPS        
---- VAR u                  0.0936         0.0965         0.0967          .         
---- VAR v                 -INF            0.0500        +INF             .         
---- VAR w                  0.0002         0.0003         0.0003         EPS        
---- VAR s                 -INF            5.3000        +INF             .         
---- VAR b1                 1.0000         1.0000         1.0000         EPS        
---- VAR b2                  .              .              .             EPS        

Tuesday, March 24, 2009

Endogenous if-then-else

> I am trying to define a nonlinear function where the variable in the
> objective is dependent on a parameter (x) which varies across
> i=1,2,...,N cases and two variables, a cut-off value (cutoff) and a
> coefficient (coeff) which transforms x above the value of cutoff. The
> logic, in pseudocode is:
>
> if x[i] > cutoff {
>     x2[i] = x[i] * coeff
>     }
> else {
>     x2[i] = x[i]
>     }
>
> The objective function would be something like:
>
> minimize Objective: sum {i in N} [complicated function goes here] * x2[i] ;
>
> How, if at all, can this be implemented in AMPL?
>

You could try:
e{i in I}: x2[i] = if (x[i]>cutoff) then x[i]*coeff else x[i];
However this is a discontinuous function so beware that many NLP solvers assume smooth functions and may get into trouble near the jump. A smooth approximation would be better in that case (there are several possible functional forms for such a smooth approximation; one would need to invest some time to try a few of them; an example is shown here). If the model was otherwise linear, one could look into a piecewise linear formulation (AMPL has very good support for this; this approach is less interesting if the model is otherwise nonlinear as it would become an MINLP).

Monday, March 23, 2009

Indexing by a variable

> I've:
>
> param N_MAX = 1000;
> var N >= 1, <= N_MAX;
>
> var x {0 .. N_MAX} >= 0;
>
> I want to do *something like*:
>
> subject to {i in  {(N/2)-1, (N/2)+1} }: #calculation of indices
> simplified for sake of discussion
>     x[i] = 1000; #RHS simplified for sake of discussion
>
> But of course I can't use variables for indexing. Neither can I do it
> explicitly as there are too many possible values that N can take. So,
> I'm looking for suggestions.
>
> Sincerely,

param N_MAX := 10;

set I := {1..N_MAX};
set J := {0..N_MAX};
set IJ{i in I} := setof{j in J: ((j >= i/2-1) and (j <= i/2+1))}(j);
display IJ;

var N >= 1, <= N_MAX;
param xup := 10000;
var x {J} >= 0, <= xup;
var b {I} binary;

minimize obj: sum{j in J} x[j];
s.t. e1: sum{i in I} i*b[i] = N;
s.t. e2: sum{i in I} b[i] = 1;
s.t. e3{i in I, j in IJ[i]}:  x[j] <= 1000*b[i] + xup*(1-b[i]);
s.t. e4{i in I, j in IJ[i]}:  x[j] >= 1000*b[i];

solve;
display N,b,x;
end;

Usually it is best to introduce a binary variable b[i] with b[i]=1 <=> i=N. Then you can introduce constraints based on b[i]. xup is here an upper bound on x[j]. Choose this number as tight as possible: it is a big-M constant. As usual it is a good idea to move as much as possible complexity away from the constraints. In this case we use an intermediate set IJ to handle a more difficult concept. Such a set can be tested and debugged before a complete model is ready to run.

MSF Examples

There were some requests to make my MS Solver Foundation examples available. I have placed them here. The examples include:

  1. TSP in OML (post)
  2. Tiling squares in OML (post1 and post2)
  3. Sudoku (post)
  4. Maxflow (post)
  5. Langford’s problem (post)
  6. Magic Squares (post)

Thursday, March 19, 2009

MIP model

Hi,
I am dealing with an optimization task as follows:
1) Tom is given a large box which contains 1000 bags of marbles.
2) Inside each bag, there are between 1 and 50 marbles.
3) Within a given bag, there are no duplicate colors of marbles, however there are duplicate colors of marbles across the bags.
4) The bags are labeled and Tom has a list of what color marbles each bag contains.
Objective:  Tom wants to collect as many UNIQUE colors of marbles as possible, but is only allowed to pick out 100 bags from the box.  Tom derives absolutely no value from getting duplicate colors of marbles.  Tell Tom which 100 bags to pick.
My question:
First of all, is this even a problem that's appropriate for linear programming? I was able to solve this quite simply on a small scale in Excel's Solver, but when I moved over to GLPK I couldn't figure out how to avoid making my objective statement non-linear.  Initially I wrote out the problem like below.  However, glpk complained about having a variable in my if/else statement in the objective function.  I'm still rather new to LP so if someone could tell me whether or not I should even be using LP for this, and if so, I'd appreciate some suggestions on how to model the problem such that GLPK doesn't complain. Thanks.
/* sets */
set COLORS;
set BAGS;
/* parameters */
param colors_to_bags {i in COLORS, j in BAGS}
/* decision variables */
var selected_bags {j in BAGS} binary >= 0;
/* objective function */
maximize uniques: sum{i in COLORS} if (sum{j in BAGS} colors_to_bags[i, j] * selected_bags[j]) > 1 then 1 else 0;
/* constraints */
s.t. max_bags: sum{j in BAGS} selected_bags{j} <= 100;
data;
set BAGS := bag1 bag2 .... bag1000;
set COLORS := green yellow blue ...etc; (assume there are many many colors)
/* 1 means the bag contains that color marble, 0 means it does not */
praram colors_to_bags: bag1 bag2 .... bag1000:=
green 0 1 0 0 ...
yellow 0 1 0 1 ...
blue 1 0 0 1 ...
etc. 0 1 1 0 ...;
end;

You need to introduce an extra variable have_color[i]:

var have_color{i in COLORS} binary;

/* all colors_to_bags[i,j]*selected_bags[j] = 0 =>have_color[i]=0 */
eq1{i in COLORS}: have_color[i] <= sum{j in BAGS} colors_to_bags[i,j]*selected_bags[j];

/* any colors_to_bags[i,j]*selected_bags[j] = 1 =>have_color[i]=1 */
eq2{i in COLORS, j in BAGS}: have_color[i] >= colors_to_bags[i,j]*selected_bags[j];

Then the obj becomes:

maximize uniques: sum{i in COLORS} have_color[i];

You can improve performance by relaxing have_color to be a continuous variable between 0 and 1 (it is automatically integer). In addition you may drop eq2 as you are maximizing this anyway. We are lucky: eq2 is adding many equations.