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
CodePudding user response：
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]
CodePudding user response：
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.