PYOMO indexed variable with different domains

[Marc Meketon]

I just "discovered" a feature of PYOMO that isn't in every modeling language, and yet is quite useful and simple.  I'm hoping by posting this note it will allow others in the future to search for this trick.

The overall concept is that a variable indexed over a set cannot have a "rule" that determines the domain for each set/variable element.  This is in contrast to the bounds of a variable, where you can have a rule that sets the bounds for each set/variable.  It might be desirable to have different domains for different indices of a variable set, especially when experimenting with integerization heuristics that force as binary only a subset of the variables that would ordinarily be binary.

Suppose that you have an model M, a set M.K, and a variable M.X indexed by M.K.  [I'm going to drop the "M." at this point.]  There are situations in which you might want some of the X[k]'s (k in K) to be binary, and others to be continuous, but you don't want to complicate the constraints with lots of different sums over different sets.  To do so:

• Break up K into two disjoint sets, KB and KC, where KB represents the indices for the binary variables and KC are the indices for the continuous variables.  Keep K as well.
  
• Break up variable X into two variables, XB and XC, where XB is defined over the index set KB and XC is over KC.  Don't keep the original X.
• Create a dictionary with keys in K and values that are either XB or XC, depending on whether the key is in KB or KC.
• In the constraints and objective function, use the dictionary X. 

Here's a simple example.  Consider the classical Knapsack problem, where you try to find a subset of items that have the maximum profit but fit within a weight capacity.  Typically (not always... this is a heuristic!), the items that have large profit to weight ratio will be included, the items with small profit to weight will be excluded.  Let's pretend (for the sake of this example) that the capacity is about one half of the total weight of all items.  A reasonable heuristic that reduces the number of binary variables might be to have all the items whose profit-to-weight ratio is in the middle 25% to be binary, and to have the items outside that set to be continuous.  After solving, throw away the (at most one) variable that is not 0/1.  I know this is a silly example, but it shows the use of the dictionary to select which variable is used.  An implementation is at the end.  Note (and this is the important part) that the constraints and objectives do not need separate terms for the binary and continuous variables.

Another example -- one that is the subject of a current experiment -- is where there are huge numbers of binary variables indexed over a set of (say) dimension 4.  One of the dimensions is "time", and let's say that "time" goes from 1 to 40.  The experiment is to have, say, the first 10 time periods be binary, the other, later 30 time periods to be continuous in [0.0,1.0].  The MIP is solved, then the first, say, 8 time periods are fixed, and now re-solved with periods 9 - 18 as binary and time periods 19 - 40 as continuous.  Re-solve, fix variables from 9 - 16, make variables 17 - 26 binary, re-solve, and so on. It is too early for me to recommend this path, but at least this should give an understanding of the usefulness of this approach and why the use of dictionary whose values are PYOMO variables is useful.

from __future__ import division

from pyomo.environ import *

from random import randint, seed

seed(321)

N = 50

weight = [randint(1,100) for i in range(0,N)]

profit = [randint(1,100) for i in range(0,N)]

totalWeight = sum(weight)

KnapsackCapacity = totalWeight * 0.55

profit_to_weight = [profit[i] / weight[i] for i in range(0,N)]

import numpy as np

lowerPercentile = np.percentile(profit_to_weight,37.5)

upperPercentile = np.percentile(profit_to_weight,62.5)

BinaryIdx = {i+1 for i in range(0,N) 

                   if lowerPercentile <= profit_to_weight[i] and profit_to_weight[i] <= upperPercentile}

ContinuousIdx = {i for i in range(1,N+1)} - BinaryIdx

m = ConcreteModel()

m.xb = Var(BinaryIdx, domain=Binary )

m.xc = Var(ContinuousIdx, bounds=(0.0,1.0))

m.I = Set(initialize = RangeSet(N))

X = {i : m.xb[i] if i in BinaryIdx else m.xc[i] for i in m.I}

m.KnapsackConstraint = Constraint(expr = sum(weight[i-1]*X[i] for i in m.I) <= KnapsackCapacity)

m.TotalProfit = Objective(expr = sum(profit[i-1]*X[i] for i in m.I), sense=maximize)

opt = SolverFactory("glpk")

results = opt.solve(m, tee=True)

print("Profit    weight    p2w    type value")

for im1 in range(0,N):

    i = im1 + 1

    if i in BinaryIdx:

        type = "B"

    else:

        type = "C"

    print('{0:7.3f}  {1:7.3f}{2:8.2f}      {3}  {4:4.2f}'.format(profit[im1], weight[im1], profit_to_weight[im1], type, X[i].value))

[Gabriel Hackebeil]

Marc,

Thanks for sharing this detailed example. Let me just comment that we do in fact support multiple domains within a single indexed variable, but it’s been a bit of a hidden feature up until now.

If you are on Pyomo 4.1 and 4.2, you can “activate” this feature by defining the variable domain using a rule.

model.X = Var(model.s, domain=lambda m,i: Reals)

This was recently changed on Pyomo trunk so that you no longer need to do this. That is, all indices of an indexed variable now have their own domain by default.

model = ConcreteModel()

model.s = Set(initialize=[1,2,3])

model.X = Var(model.s, domain=Binary)

model.X[2].domain = UnitInterval

This functionality is not guaranteed to be in the next release yet, but the original way of activating it by declaring the domain using a rule would probably stick around either way.

Gabe

--
You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Marc Meketon]

Thank you Gabe for pointing this out.   I notice your use of UnitInterval.  The PYOMO documentation (both the book and online documentation such as https://software.sandia.gov/downloads/pub/pyomo/PyomoOnlineDocs.html#_predefined_virtual_sets  that I have seen has not listed this as a predefined set, although PercentFraction (which seems to be the same [0,1] interval) is there. Is there more up-to-date documentation available?

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "Pyomo Forum" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/pyomo-forum/1hgQEhrjEeE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to pyomo-forum+unsubscribe@googlegroups.com.

[Gabriel Hackebeil]

The most up to date docs can be found online at pyomo.org. However, these can be both behind and in front of the latest release. We try to keep them up to date as things change on Pyomo trunk (but forget sometimes), and, as far as I can tell, this is the documentation that we link to on that site.

I added a line about UnitInterval there, so thanks for pointing that out.

Gabe

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.