Linear programming/MIP setting up sum of conditional constraints?

Pewter City 0
Viewed 39

I'm trying to set up MIP problem using LpSolve and I'm at a loss on how to set it up.

We are trying to allocate resources to various investments, to maximize revenue, subject to a budget constraint. Each business is located in a particular state, and We also have to spread our allocations across 5 states (this is the one that's causing me issues). My current linear equations are this:

Objective function: maximize following where Xi represents the number of dollars allocated to ith business, and Ri represents revenue generated by the ith business

X1 * R1 +....+ X35 * R35

Subject to the constraint:

X1 * C1 +....+ X35 * C35 <= 1,000,000

Where Ci represents the costs of the ith business. Setting this objective function and constraint matrix is very easy. However each business is located in some state, and what I need to do is setup a constraint that we must allocate to business in at least 5 states.

What I initially considered was a set of constraints where each state was represented as a row: so if this row represents the state of California for instance:

X1 * 1 +....+ X2 * 1 + X3 *0 + X35 * 1 <= Some Constraint

Then each business in California is assigned a value of 1, then this would give me the sum of $s (from the Xi assigned to each business) allocated to California. I could do this for all states, but that doesn't quite get me what I want as I thought about it mroe. I could constrain each state using this methodology, but what I want is something that lets me assign Unique_States >5

But this doesn't quite get me what I want. What I really want is something that gives me a sum of total number of states "used".

Is this even possible?

I've spent a while trying to rack my brain on how to do this but I can't quite figure it out.

There's one other questions like this on SO, but I couldn't follow any of the answers: LpSolve R conditional constraint

This question is awfully similar albeit a totally different context, but I couldn't follow it and I don't have enough reputation to ask any questions.

If anyone could give any sort of guidance

lpsolve mixed-integer-programming linear-programming constraints r
2 Answers

Here is a skeleton of your model in pyomo. This is basically a "Knapsack Problem" with a twist (the min number of states reqt.).

Code

import pyomo.environ as pyo
from collections import defaultdict

# some data...  this could be from dictionary (as below) or file read or ....

#               name                 state  C   R
investments = { 'Burger Joint':     ('CA', 1.2, 3.1),
                'car wash':         ('CA', 1.1, 2.6),
                'gem store':        ('NV', 1.0, 2.5),
                'warehouse':        ('UT', 1.2, 3.6)}

budget = 10
min_investment = 1
min_states_to_use = 2

states = {v[0] for v in investments.values()}
investments_by_state = defaultdict(list)
for k, v in investments.items():
    investments_by_state[v[0]].append(k)

# set up the model

m = pyo.ConcreteModel()

# SETS

m.I = pyo.Set(initialize=investments.keys())
m.S = pyo.Set(initialize=list(states))
# the below is an "indexed set" which is indexed by STATE and contains the 
# subset of I within that state
m.IS = pyo.Set(m.S, within=m.I, initialize=investments_by_state)

# PARAMETERS
m.cost = pyo.Param(m.I, initialize={k:v[1] for k, v in investments.items()})
m.rev  = pyo.Param(m.I, initialize={k:v[2] for k, v in investments.items()})

# VARIABLES
m.X = pyo.Var(m.I, domain=pyo.NonNegativeReals)  # amount invested in i
m.Y = pyo.Var(m.S, domain=pyo.Binary)            # 1 if state y has an investment

# OBJECTIVE
m.obj = pyo.Objective(expr=sum(m.X[i]*m.rev[i] for i in m.I), sense=pyo.maximize)

# CONSTRAINTS
# stay in budget...
m.C1 = pyo.Constraint(expr=sum(m.X[i]*m.cost[i] for i in m.I) <= budget)

# connect indicator Y to X...
def state_invest(m, s):
    return m.Y[s] * min_investment <= sum(m.X[i] for i in m.IS[s])
m.C2 = pyo.Constraint(m.S, rule=state_invest)

# use needed number of states
m.C3 = pyo.Constraint(expr=sum(m.Y[s] for s in m.S) >= min_states_to_use)

m.pprint()  # <-- show the whole model

# instantiate a solver and solve...
solver = pyo.SolverFactory('cbc')
result = solver.solve(m)
print(result)

m.X.display()  # <-- "display" substitutes in the variable values
m.Y.display()

Output

3 Set Declarations
    I : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    4 : {'Burger Joint', 'car wash', 'gem store', 'warehouse'}
    IS : Size=3, Index=S, Ordered=Insertion
        Key : Dimen : Domain : Size : Members
         CA :     1 :      I :    2 : {'Burger Joint', 'car wash'}
         NV :     1 :      I :    1 : {'gem store',}
         UT :     1 :      I :    1 : {'warehouse',}
    S : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    3 : {'CA', 'NV', 'UT'}

2 Param Declarations
    cost : Size=4, Index=I, Domain=Any, Default=None, Mutable=False
        Key          : Value
        Burger Joint :   1.2
            car wash :   1.1
           gem store :   1.0
           warehouse :   1.2
    rev : Size=4, Index=I, Domain=Any, Default=None, Mutable=False
        Key          : Value
        Burger Joint :   3.1
            car wash :   2.6
           gem store :   2.5
           warehouse :   3.6

2 Var Declarations
    X : Size=4, Index=I
        Key          : Lower : Value : Upper : Fixed : Stale : Domain
        Burger Joint :     0 :  None :  None : False :  True : NonNegativeReals
            car wash :     0 :  None :  None : False :  True : NonNegativeReals
           gem store :     0 :  None :  None : False :  True : NonNegativeReals
           warehouse :     0 :  None :  None : False :  True : NonNegativeReals
    Y : Size=3, Index=S
        Key : Lower : Value : Upper : Fixed : Stale : Domain
         CA :     0 :  None :     1 : False :  True : Binary
         NV :     0 :  None :     1 : False :  True : Binary
         UT :     0 :  None :     1 : False :  True : Binary

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : 3.1*X[Burger Joint] + 2.6*X[car wash] + 2.5*X[gem store] + 3.6*X[warehouse]

3 Constraint Declarations
    C1 : Size=1, Index=None, Active=True
        Key  : Lower : Body                                                                    : Upper : Active
        None :  -Inf : 1.2*X[Burger Joint] + 1.1*X[car wash] + X[gem store] + 1.2*X[warehouse] :  10.0 :   True
    C2 : Size=3, Index=S, Active=True
        Key : Lower : Body                                    : Upper : Active
         CA :  -Inf : Y[CA] - (X[Burger Joint] + X[car wash]) :   0.0 :   True
         NV :  -Inf :                    Y[NV] - X[gem store] :   0.0 :   True
         UT :  -Inf :                    Y[UT] - X[warehouse] :   0.0 :   True
    C3 : Size=1, Index=None, Active=True
        Key  : Lower : Body                  : Upper : Active
        None :   2.0 : Y[CA] + Y[NV] + Y[UT] :  +Inf :   True

11 Declarations: I S IS cost rev X Y obj C1 C2 C3

Problem: 
- Name: unknown
  Lower bound: -29.5
  Upper bound: -29.5
  Number of objectives: 1
  Number of constraints: 5
  Number of variables: 7
  Number of binary variables: 3
  Number of integer variables: 3
  Number of nonzeros: 4
  Sense: maximize
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.0
  Wallclock time: 0.0
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
    Black box: 
      Number of iterations: 0
  Error rc: 0
  Time: 0.014362096786499023
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

X : Size=4, Index=I
    Key          : Lower : Value     : Upper : Fixed : Stale : Domain
    Burger Joint :     0 :       1.0 :  None : False : False : NonNegativeReals
        car wash :     0 :       0.0 :  None : False : False : NonNegativeReals
       gem store :     0 :       0.0 :  None : False : False : NonNegativeReals
       warehouse :     0 : 7.3333333 :  None : False : False : NonNegativeReals
Y : Size=3, Index=S
    Key : Lower : Value : Upper : Fixed : Stale : Domain
     CA :     0 :   1.0 :     1 : False : False : Binary
     NV :     0 :   0.0 :     1 : False : False : Binary
     UT :     0 :   1.0 :     1 : False : False : Binary
[Finished in 696ms]
AirSquid0

There are 2 general approaches here that might help. Caveat: I don't know LpSolve syntax and I don't know how you would (or could) express these in r/LpSolve syntax. Aside: Are you open to a different framework? There are several LP frameworks in python that would make this very easy to set up... Anyhow... with some pseudocode...

In both of the approaches below, you will need to implement a binary "indicator variable" that indicates whether a state has been used. This also assumes there is some logical "min investment" value to prevent spreading something like 0.00001 in a particular investment

S = {CA, NV, UT, ...}  set of states
Y[s ∈ S] ∈ {0, 1}      binary indicator

Approach 1: Make subsets of your investment options by state and index those by state. Then you could sum over the whole set as needed

I = {1, 2, ..., 35}              set of investment options
X[i ∈ I] ∈ {non-negative reals}  the amount invested in i
IS = {CA: {1, 5, 10,...}, NV: {2, 4, 15,...} , ...}   subsets by state

for each s ∈ S:
    Y[s] * min_investment <= ∑ X[IS[s]]   # constrain Y to be less than state sum

∑Y[s] >= 5

Approach 2: double-index your investments by [state, i] and then you can sum over either the states or investments (I).

X[s, i] <= 0   constrain for all applicable combos that are illegal
for s ∈ S:
    Y[s] * min_investment <= ∑ (X[s, i] for all i ∈ I)
∑Y[s] >= 5

It appears that LpSolve wants a big matrix for constraints, so either of these might be difficult in that framework, but they are straightforward in other frameworks. Perhaps somebody w/ deeper experience in LpSolve can suggest syntax.

AirSquid0
Related