1

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

Improve this question

2 Answers 2

Reset to default
2

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]
Improve this answer
Sign up to request clarification or add additional context in comments.

1 Comment

Can't thank you enough. This has been bugging me for 10 days now, thank you very very very much, I really appreciate your generosity in the answers. Above and beyond anything I could have asked for. Thank you!
0

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.

Improve this answer

3 Comments

More than happy to investigate other models other than lpSolve. I looked into glpk but it seemed to have awfully similar syntax to lpSolve. Your pseudocode that begins with "for each s in S" is absolutely perfect, I had no idea you could use syntax like that in a linear optimizer since I've only ever used lpSolve which only allows numerical matrix inputs (from what I can tell)
Do you have packages/libraries you could recommend? Googling suggests PULP might be a popular one, is that one that would work similar to the pseudocode you laid out above?
I use both pulp and pyomo. I prefer pyomo. It is more expressive and the syntax is clearer, IMHO. Both use fairly straightforward python syntax to build optimization models. Pyomo does NOT include a solver, but there are several free ones you can install and use any with pyomo. GLPK, CBC are both good (free) solvers for LP. I've posted quite a few answers in both pyomo and pulp tags.... take a look at some of my answers for ideas. I'll code you an example for your problem above in pyomo......stby...

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.