In my optimization problem I have to allocate as much wagons as possible to its optimal routes for minimum cost with defined constraints (multi objective problem). At first step I have to distribute as much wagons as possible (maximum amount) and the second step of optimization is to allocate its number to the cheapest routes. I have the number of wagons from each station to distribute and I have the capacity at each station of arrival. For simplicity, I have only two stations of departure where rem_type column means the type of wagon at this station. Table df_vrp contains information about the capacity and table df_unload contains information about the number of wagons at each departure station which we have to allocate to routes.It is clear for me how to the maximum possible number of wagons from station to station in regards with constraints and then find the cheapest routes. But, in my case, I have to take in account the proportion of distributes wagons such that the sum of distriubuted wagons of Type 2 must be equal to 1/3 of the sum of distributed wagons of Type 1. So I cannot understand how to create such constraint because 1/3 can be done in many ways (1 wagon of Type 2 and 3 wagons of Type 1, 3 wagons of Type 2 and 9 wagons of Type 1 and so on).
import pandas as pd
import numpy as np
from pyomo.environ import *
data = {'rem_type' : ['Type1', 'Type1', 'Type2', 'Type2', 'Type1', 'Type1', 'Type2', 'Type2'],
'station_depart': ['A', 'A', 'A', 'A', 'F', 'F', 'F', 'F'],
'station_arrive': ['B', 'C', 'B', 'C', 'B', 'C', 'B', 'C'],
'cost': [100, 103, 111, 101, 105, 114, 95, 99]}
df_route = pd.DataFrame(data)
# unload
df_unload = pd.DataFrame({'rem_type' : ['Type1', 'Type2', 'Type1', 'Type2'],
'station_depart' : ['A', 'A', 'F', 'F'],
'volume' : [5, 6, 4, 7]})
# repair
df_vrp = pd.DataFrame({'rem_type' : ['Type1', 'Type2'],
'station_arrive' : ['B', 'C'],
'capacity' : [15, 30]})
# create routes
routes_unload = dict()
routes_vrp = dict()
# create model
model = ConcreteModel("OP")
# create multi objective function
model.Models = ObjectiveList()
# decision vars
model.x = Var([idx for idx in df_route.index], domain=NonNegativeIntegers)
# obj function to maximize the number of allocated wags to routes
model.Size = Objective(expr=sum([model.x[i] for i in model.x]), sense=maximize)
# obj function to minimize cost of those wags from the first solution
model.Cost = Objective(
expr=sum([df_route.loc[idx, "cost"] * model.x[idx] for idx in df_route.index]),
sense=minimize,
)
model.Size.activate()
model.Cost.deactivate()
# routes with indexes for creating constraints
for idx in tqdm(df_route.index):
vrp = (df_route.loc[idx, "rem_type"], df_route.loc[idx, "station_arrive"])
if vrp not in routes_vrp:
routes_vrp[vrp] = list()
routes_vrp[vrp].append(idx)
t_unload = (df_route.loc[idx, "rem_type"], df_route.loc[idx, "station_depart"])
if t_unload not in routes_unload:
routes_unload[t_unload] = list()
routes_unload[t_unload].append(idx)
# constraints on the arrive/repair station
model.vrp = ConstraintList()
for v in df_vrp.index:
t = (df_vrp.loc[v, "rem_type"], df_vrp.loc[v, "station_arrive"])
if t in routes_vrp:
model.vrp.add(
sum([model.x[idx] for idx in routes_vrp[t]]) <= df_vrp.loc[v, "capacity"]
)
# constraints on the depart station (how many wagons we have for allocation)
model.unload = ConstraintList()
for u in df_unload.index:
t = (df_unload.loc[u, "rem_type"], df_unload.loc[u, "station_depart"])
if t in routes_unload:
model.unload.add(
sum([model.x[idx] for idx in routes_unload[t]]) <= df_unload.loc[u, "volume"]
)
results = SolverFactory("glpk").solve(model)
results.write()
# the result from the first objective function as the constraint to the second one
model.fix_count = Constraint(expr=sum([model.x[i] for i in model.x]) == model.Size())
model.Size.deactivate()
model.Cost.activate()
results = SolverFactory("glpk").solve(model)
results.write()
# create df with results of allocation to the routes
vals = []
for i in tqdm(model.x):
vals.append(model.x[i].value)
df_results = df_route.dropna(
subset=["rem_type", "station_depart", "station_arrive", "cost"],
how="all",
)
df_results["total"] = vals
Here is a complete answer, based on your updated code. I went a little "geek" on this, but in thinking it over, you helped me with a different problem I've been working on, so we're GTG.
I think one of the main issues with your base code is that you are trying to force pandas into a pyomo box and it gets way difficult very quick and the code becomes pretty confusing. I'd recommend ditching pandas as quickly as possible for clarity and your ability to work with sets more cleanly because as you were doing it before, you have were indexing by the enumeration of all the routes, and you'd really like to be able to make constraints based on things like destinations, types of vehicles, etc. which is wayyy hard to do if you don't set it up that way.
I also brought all of the parameters into the model. This is optional, but I think a good practice to be able to pprint() the whole model or parts of it more cleanly.
If you have a non-complete set of pairs of routes and types, there is a bit more work to do. Meaning, if there are routes that only host 1 type of vehicle, you'd have to do a little filtering or use indexed sets or such, but right now all types are compatible on all routes, so this works.
Note the result of the 2nd pass flips some of the results to minimize cost, which is your target.
import pandas as pd # <--- try to avoid this, unless needed to curate the data
# import numpy as np
from pyomo.environ import *
data = {'rem_type' : ['Type1', 'Type1', 'Type2', 'Type2', 'Type1', 'Type1', 'Type2', 'Type2'],
'station_depart': ['A', 'A', 'A', 'A', 'F', 'F', 'F', 'F'],
'station_arrive': ['B', 'C', 'B', 'C', 'B', 'C', 'B', 'C'],
'cost': [100, 103, 111, 101, 105, 114, 95, 99]}
df_route = pd.DataFrame(data)
costs = df_route.set_index(['station_depart', 'station_arrive', 'rem_type']).to_dict()['cost']
# unload
df_unload = pd.DataFrame({'rem_type' : ['Type1', 'Type2', 'Type1', 'Type2'],
'station_depart' : ['A', 'A', 'F', 'F'],
'volume' : [5, 6, 4, 7]})
demands = df_unload.set_index(['station_depart', 'rem_type']).to_dict()['volume']
# repair
df_vrp = pd.DataFrame({'rem_type' : ['Type1', 'Type2'],
'station_arrive' : ['B', 'C'],
'capacity' : [15, 30]})
capacity = df_vrp.set_index(['station_arrive', 'rem_type']).to_dict()['capacity']
# find the set of available routes and others conveniences
routes = {(s, t) for (s, t, cost) in costs.keys()}
start_nodes = {k[0] for k in costs.keys()}
end_nodes = {k[1] for k in costs.keys()}
all_nodes = start_nodes | end_nodes
# Target Ratios
ratios = {('Type1', 'Type2') : (3, 1)
# others
}
# MODEL
# create model
model = ConcreteModel("OP")
# SETS
model.N = Set(initialize=sorted(all_nodes))
model.S = Set(within=model.N, initialize=sorted(start_nodes)) # the set of start nodes
model.T = Set(within=model.N, initialize=sorted(end_nodes)) # the set of terminations
model.R = Set(within=model.N * model.N, initialize=sorted(routes)) # the set of connections
model.V = Set(initialize=['Type1', 'Type2']) # vehicle types
model.ratio_pairs = Set(within=model.V * model.V, initialize=ratios.keys())
# PARAMS
model.cost = Param(model.R, model.V, initialize=costs)
model.demand = Param(model.S, model.V, initialize=demands, default=0) # the demand to send from s of type v
model.capacity = Param(model.T, model.V, initialize=capacity, default=0) # the capacity to accept at t type v
model.ratio_limits = Param(model.ratio_pairs, initialize=ratios, domain=Any)
# VARS
model.send = Var(model.R, model.V, domain=NonNegativeIntegers) # send this number of vehicles of type v on route r
# OBJ1: Route as much as possible
model.obj1 = Objective(expr=sum_product(model.send), sense=maximize)
# CONSTRAINTS
@model.Constraint(model.S, model.V)
def demand_limit(model, s, v):
return sum(model.send[r, v] for r in model.R if r[0] == s) <= model.demand[s, v]
@model.Constraint(model.T, model.V)
def capacity_limit(model, t, v):
return sum(model.send[r, v] for r in model.R if r[1] == t) <= model.capacity[t, v]
# the ratio constraint...
@model.Constraint(model.ratio_pairs)
def ratio_limit(model, v1, v2):
return sum(model.send[r, v1] for r in model.R) * model.ratio_limits[v1, v2][0] == \
sum(model.send[r, v2] for r in model.R) * model.ratio_limits[v1, v2][1]
# solve first pass...
results = SolverFactory("glpk").solve(model)
results.write()
model.send.display()
model.ratio_limit.pprint()
# now, add a constraint to make sure we deliver at least as much as before and introduce a
# new cost objective
model.min_send = Constraint(expr=sum_product(model.send) >= value(model.obj1))
model.obj2 = Objective(expr=sum_product(model.cost, model.send), sense=minimize)
model.obj1.deactivate()
results = SolverFactory("glpk").solve(model)
results.write()
model.send.display()
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: 16.0
Upper bound: 16.0
Number of objectives: 1
Number of constraints: 10
Number of variables: 9
Number of nonzeros: 25
Sense: maximize
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Termination condition: optimal
Statistics:
Branch and bound:
Number of bounded subproblems: 9
Number of created subproblems: 9
Error rc: 0
Time: 0.006043910980224609
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
send : Size=8, Index=send_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
('A', 'B', 'Type1') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('A', 'B', 'Type2') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('A', 'C', 'Type1') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('A', 'C', 'Type2') : 0 : 5.0 : None : False : False : NonNegativeIntegers
('F', 'B', 'Type1') : 0 : 4.0 : None : False : False : NonNegativeIntegers
('F', 'B', 'Type2') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('F', 'C', 'Type1') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('F', 'C', 'Type2') : 0 : 7.0 : None : False : False : NonNegativeIntegers
ratio_limit : Size=1, Index=ratio_pairs, Active=True
Key : Lower : Body : Upper : Active
('Type1', 'Type2') : 0.0 : (send[A,B,Type1] + send[A,C,Type1] + send[F,B,Type1] + send[F,C,Type1])*3 - (send[A,B,Type2] + send[A,C,Type2] + send[F,B,Type2] + send[F,C,Type2]) : 0.0 : True
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: 1598.0
Upper bound: 1598.0
Number of objectives: 1
Number of constraints: 11
Number of variables: 9
Number of nonzeros: 33
Sense: minimize
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Termination condition: optimal
Statistics:
Branch and bound:
Number of bounded subproblems: 1
Number of created subproblems: 1
Error rc: 0
Time: 0.00722813606262207
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
send : Size=8, Index=send_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
('A', 'B', 'Type1') : 0 : 4.0 : None : False : False : NonNegativeIntegers
('A', 'B', 'Type2') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('A', 'C', 'Type1') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('A', 'C', 'Type2') : 0 : 5.0 : None : False : False : NonNegativeIntegers
('F', 'B', 'Type1') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('F', 'B', 'Type2') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('F', 'C', 'Type1') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('F', 'C', 'Type2') : 0 : 7.0 : None : False : False : NonNegativeIntegers
[Finished in 631ms]