Vectorized constraint

回答済み

Bushido Pora

2023年04月20日 11:09

• 
• 
• 

Hello,

Im trying to encode a constraint in a vectorized way because for loops are far too slow

Here is the classic constraint :

for i in range(10):

    for i_tilde in range(10):

        for j in range(20):

            model.addConstr(x[i,j]+x[i_tilde,j] <= 1)

Since the constraint stands for all j, it seems natural to vectorize it as A * (x[i,:]+x[i_tilde,:]) <= b with A = np.ones((20,20)) and b = np.ones(20)

However I could not find a way to code this in gurobi, could you help me out ?

I'm looking for something like :

model.addMConstr(A = np.identity(20, dtype=int), x = x[i,:]+x[i_tilde,:], sense = '<=', b = np.ones(20) )

Thanks!

0

• 

  Silke Horn

  • Gurobi Staff

  2023年04月20日 13:06

  Hi,

  For me, the most natural way to implement those constraints would actually be the addConstrs method:

  x = model.addVars(10, 20, name="x")
  model.addConstrs(x[i,j] + x[i_tilde,j] <= 1 for i in range(10) for i_tilde in range(10) for j in range(20))

  Note, though, that this does not yield any performance benefits over the for-loops. However, with these dimensions, it does not take any noticeable time on my laptop. How big is your model going to be?

  I don't see any natural way to implement this using matrices.

  - Silke

  0
• 

  Riley Clement

  • Gurobi Staff

  2023年04月21日 00:47

  Hi Bushido,

  Just to confirm, your x variables are continuous, not binary, right?

  - Riley

  0
• 

  Bushido Pora

  2023年04月24日 07:28

  • 
  • 
  • 

  Hi,

  No, x is binary and has quite a large size (20000,500) thats why the for loop is huge, it performs 20k²*500 = 200000000000 iterations which is intractable. 

  That is why I need a matrix form, 

  Thanks in advance

  0
• 

  Riley Clement

  • Gurobi Staff

  2023年04月24日 08:29

  Hi Bushido,

  If x is binary you will need to tweak what you are doing to ensure correctness.  When i == i_tilde you will have

  x[i,j]+x[i,j] <= 1

  for all j, which will necessarily mean that x[i,j] cannot ever be 1 and all your x variables will be zero.

  You can rectify this by making sure you only generate the constraints for i != i_tilde.  However there is an opportunity to make your formulation stronger, and the model building faster.

  If you have the following three constraints for some j

  x[0,j]+x[1,j] <= 1
  x[0,j]+x[2,j] <= 1
  x[1,j]+x[2,j] <= 1

  then only one of these x variables can have a value of one (test it for yourself to understand why).  This means we can replace these 3 constraints with what is known as a clique constraint:

  x[0,j]+x[1,j]+x[2,j] <= 1

  Likewise you can replace your constraints with

  model.addConstrs(gp.quicksum(x[i,j] for i in range(10)) <= 1 for j in range(20))

  or with matrix variables (note j is handled automatically)

  model.addConstr(x.sum(axis=0) <= 1, name="con_j=")

  - Riley

  0
• 

  Bushido Pora

  2023年04月24日 12:57 編集済み

  • 
  • 
  • 

  Hello,

  My bad for not giving more contexts, I'm modelling a problem similar to the aircraft allocation problem :

  x(i,j) = {1 if flight i takes place in aircraft j, 0 otherwise}

  let INC(i) denote the indexes of flights which have a non-empty schedule intersection with flight i (from which we remove i itself)

  Since any flight in INC(i) cannot use the same aircraft as i: x_{i,j} + x{i_tilde, _j} <=1 for all i_tilde in INC(i)

  Encoding this constraint seemed to require a 3-nested loop to me.

  Thanks for answer above.

  Maybe I misunderstood your example but in this setting one cannot use a summation constraint as 2 compatible flights might be cancelled together because each would be incompatible with a third one.

  Suppose flight 12 is incompatible with flight 11 and flight 10 but the 2 latter are compatible together the summation constraint states that x[12,j]+x[11,j]+x[10,j] <= 1 and thus only one flight can use aircraft j even if 

  11 and 10 were compatible. The formulation pair by pair prevents those undesired behaviours.

  Note I might have missed your point.

  Thanks for your help.

  0
• 

  Riley Clement

  • Gurobi Staff

  2023年04月24日 14:15

  Hi Bushido,

  More generally, you can think of a graph where there is a node for each variable x[i,j], and if two variables are incompatible then an edge exists between them in the graph. In this graph you can find cliques, which are sets of nodes which induce complete graphs, that is, every node is connected to every other by an edge.  If we have a clique then you can form a clique constraint by restricting the sum of corresponding x variables to be less than or equal to 1.  This is a much stronger constraint than a family of constraints for each edge in the clique.

  This concept can still be applied to your problem, if there is logic that makes it easy to identify groups of flights/aircraft allocations which are all pairwise incompatible.  Although it sounds like speed will be a concern even if the identifying cliques was easy. 

  In your original example looping over range(10) twice for the flights and forming the constraints implied that flights 0..9 are all incompatible with each other, and form a clique, but also you would have had flights being incompatible with themselves.  Hopefully i never belongs to INC(i) in your definitions.

  In terms of speed of model building, as Silke noted there is no natural way of building the constraints with mathematical matrix manipulations, however we can appeal to numpy indexing tricks.  Consider the following example to motivate the idea

  >> import numpy as np
  >> arr = np.reshape(range(12), newshape=(4,3))
  >> arr
  array([[ 0,  1,  2],
         [ 3,  4,  5],
         [ 6,  7,  8],
         [ 9, 10, 11]])
  
  >> arr[[2,0,2,1,3]]
  array([[ 6,  7,  8],
         [ 0,  1,  2],
         [ 6,  7,  8],
         [ 3,  4,  5],
         [ 9, 10, 11]])

   

  Now say you have two list-like structures, eg lists, tuples or columns of a numpy array, and when paired they give the incompatibilities.  For example if you had a list INC of (i,i_tilde) pairs for incompatibilities then you can create these two iterables with the following:

  it1 = [i for i,_ in INC]
  it2 = [i for _,i in INC]

  You can then create all incompatibility constraints (for all j) at once with

  m.addConstr(x[it1] + x[it2] <= 1)

  The key to speed has now moved to how fast you can create it1 and it2, which will depend on how the incompatibility data is stored.  The itertools module, in particular chain functions, may help.  Note that you should try to avoid symmetries with incompatibilities, eg you don't need both

  x[0,3] + x[2,3] <= 1
  x[2,3] + x[0,3] <= 1

  as they are the same constraint.

  - Riley

  0

サインインしてコメントを残してください。