Passing exact number of allocations as constraint to pyomo in a sourcing problem

Question

I am solving a sourcing allocation optimization problem. Here I have let's say two brands. Each brand has a raw material demand across the 3 plants (Demand in kg)

	Brand 1	Brand2
Plant 1	3000	2000
Plant 2	2000	1000
Plant 3	3000	2000

This demand has to be satisfied by allocating vendors to fulfill demand s.t. minimize the cost of sourcing from them.

• Vendor 1 : 5000 (capacity in kg)
• Vendor 2 : 5000
• Vendor 3 : 5000
• Vendor 4 : 5000

I also have a cost matrix that suggests

Brand1 : plant 1 : vendor 1 : 20 (cost per kg)

I am struggling to define a model with a constraint as follows:

1. Allocate exactly two vendors for every brand

This is of course along with constraints

1. Non-negative quantity
2. Total allocated quantity meets the demand requirement
3. Total allocated quantity is less than equal to the capacity of the vendor

Here is the piece of code that I have tried. I am not able to create the equation correctly in pyomo which ensures a total number of vendors allocated at brand level is == 2. For the same brand, it can be one for a plant or repeated across plants.

demand = {('B1','P1'):5000,('B1','P2'):3000,('B2','P1'):5000,('B2','P2'):3000}
vendor = {'v1':2000,'v2':10000,'v3':6000}
plant = {'P1':6000,'P2':10000}
brand = {'B1':8000,'B2':8000}
cost = {('B1','v1','P1'):19,('B1','v1','P2'):27,('B1','v2','P1'):23,('B1','v2','P2'):27,('B1','v3','P1'):20,('B1','v3','P2'):20,
        ('B2','v1','P1'):19,('B2','v1','P2'):27,('B2','v2','P1'):23,('B2','v2','P2'):27,('B2','v3','P1'):20,('B2','v3','P2'):20}
model = ConcreteModel()
model.dual = Suffix(direction=Suffix.IMPORT)

model.i = Set(initialize=list(brand.keys()), doc='Brand')
model.j = Set(initialize=list(vendor.keys()), doc='Vendors')
model.k = Set(initialize=list(plant.keys()), doc='Plant')

model.a = Param(model.j, initialize=vendor, doc='Capacity of supplier i in cases')
model.b = Param(model.i,model.k, initialize=demand, doc='Demand at plant (B,P) (i,k) in cases')


model.c = Param(model.i, model.j,model.k, initialize=cost, doc='Transport cost in thousands of dollar per case')

model.x = Var(model.i, model.j,model.k, bounds=(0.0,None), doc='Shipment quantities in case')

model.y = Var(model.i, model.j,model.k, bounds=(0.0,1.0), doc='max allocation per brand')

def supply_rule(model, j):
  return sum(model.x[i,j,k] for i in model.i for k in model.k) <= model.a[j]
model.supply = Constraint(model.j, rule=supply_rule, doc='Observe supply limit at plant i')

def demand_rule(model, i,k):
  return sum(model.x[i,j,k] for j in model.j) == model.b[i,k]  
model.demand = Constraint(model.i,model.k, rule=demand_rule, doc='Satisfy demand at market j')

// ** This is the constraint i am struggling to code **//

model.count_con6 = ConstraintList()
for i in model.i:
  for j in model.j:
    if model.j!=model.j[-1]: 
      model.count_con6.add(sum(model.y[i,j,k] for k in model.k)==2)

def objective_rule(model):
  return sum(model.c[i,j,k]*model.x[i,j,k]**model.y[i,j,k] for i in model.i for j in model.j for k in model.k)
model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')

//** The above model definition returns quadratic equation error in cbc solver and model errr in ipopt**//

if __name__ == '__main__':
    # This emulates what the pyomo command-line tools does
    from pyomo.opt import SolverFactory
    import pyomo.environ
    #opt = SolverFactory("mindtpy")
    #results = opt.solve(model)
    #sends results to stdout
    SolverFactory('ipopt', executable='/content/ipopt').solve(model)
    results.write()
    print("\nDisplaying Solution\n" + '-'*60)
    pyomo_postprocess(None, model, results)

Answer 1

I will first formalize the described problem in mathematical notations, I think it is more suitable, rather than a piece of code, for providing a better understanding.

So, let

• $B$ be the brands set;
• $P_b$ be the set of plants belonging to the brand $b \in B$;
• $d_p \in \mathbb{R}^{+}$ be the plant $p \in P_b$ demand, such that $b \in B$;
• $V$ be the vendors set;
• $c_v \in \mathbb{R}^{+}$ be the vendor $v \in V$ capacity; and
• $s_{vpb} \in \mathbb{R}^{+}$ be the vendor $v \in V$ cost to supply the plant $p \in P_b$ belonging to the brand $b \in B$.

Note that, since you have not specified whether the input magnitudes are discrete, I used real domains for the parameters, however, in case you are working with integral values you can replace $\mathbb{R}^{+}$ by $\mathbb{N}$.

Now, let's go to our mathematical programming model. The variables and the two first constraints of this model are exactly equal to the model presented in your code snippet. The differences between my model and yours, are in the objective function and in third and so on constraints. So, concerning the variables, let

• $x_{vpb} \in \mathbb{R}^{+}_{\leqslant M_{vp}}$ be the amount of supply from vendor $v \in V$ allocated to plant $p \in P_b$, such that $b \in B$; and
• $y_{vb} \in \mathbb{B}$ be the flags that tell whether vendor $v \in V$ was allocated to brand $b \in B$.

In this $x$ variable, we have the big-M $M_{vp} = \text{min}\{c_v, d_p\}$ representing the maximum of supply that vendor $v \in V$ can deliver to plant $p \in \bigcup_{b \in B} P_b$. This constant also will be used in a further constraint.

The objective function.

$\text{min }\sum_{v \in V} \sum_{b \in B} \sum_{p \in P_b} x_{vpb} s_{vpb}$

Stating that each plant has its demand satisfied (demand_rule).

$\sum_{v \in V} x_{vpb} = d_p$ $\forall b \in B, p \in P_b$

Respecting the vendors' capacities (supply_rule).

$\sum_{b \in B} \sum_{p \in P_b} x_{vpb} \leqslant c_v$ $\forall v \in V$

From now on we have the required constraints you have not added to your model. Updating $y_{vb}$ whenever vendor $v \in V$ is allocated to some plant $p \in P_b$ of $b \in B$.

$y_{vb} M_{vp} \geqslant x_{vpb}$ $\forall v \in V, b \in B, p \in P_b$

Also, we have to update $y_{vb}$ whenever vendor $v \in V$ is not allocated to any plant $p \in P_b$ of $b \in B$. In this constraint, $\epsilon$ represents a minimum of demand that vendor $v$ must deliver if it is allocated to a brand $b \in B$. You can use your solver degree of tolerance value, or just some small floating number, e.g., 1e-4.

$y_{vb} \epsilon \leqslant \sum_{p \in P_b} x_{vpb}$ $\forall v \in V, b \in B$.

Case you are working with $x_{vpb} \in \mathbb{N}_{\leqslant M_{vp}}$ $\forall v \in V, b \in B, p \in P_b$, then you can remove the $\epsilon$-term, leading to the inequality below.

$y_{vb} \leqslant \sum_{p \in P_b} x_{vpb}$ $\forall v \in V, b \in B$

And finally, we force every brand to take exactly two vendors.

$\sum_{v \in V} y_{vb} = 2$ $\forall b \in B$

If anyone has any doubts or questions, let me know.