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,
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)
    print("\nDisplaying Solution\n" + '-'*60)
    pyomo_postprocess(None, model, results)
  • 1
    $\begingroup$ Welcome to OR Stack Exchange. The site supports writing mathematical notation using MathJax. Writing your model in MathJax rather than displaying Pyomo code will allow people who do not use Pyomo to answer. $\endgroup$
    – prubin
    Commented Apr 11, 2022 at 15:29
  • $\begingroup$ Are you assuming that a plant will keep track of which raw material comes from which vendor and allocate raw material to brand orders based on the vendor? $\endgroup$
    – prubin
    Commented Apr 11, 2022 at 15:30
  • $\begingroup$ thanks for the response prubin...what you described isnt a process..all the raw materials from all the vendors will be treated equally at a plant level for now...what i am struggling with is passing a count of vendor allocations at brand level as constraint restricting to exactly 2 $\endgroup$ Commented Apr 12, 2022 at 5:33
  • $\begingroup$ If a plant buys from all three vendors and treats all raw materials as equal, and that plant produces both brand 1 and brand 2, how can you know whether brand 1 is using exactly two vendors versus all three? $\endgroup$
    – prubin
    Commented Apr 12, 2022 at 18:05
  • $\begingroup$ @LakhotiaDipesh, would you elaborate more on the problem data and what you want? $\endgroup$
    – A.Omidi
    Commented Apr 14, 2022 at 19:28

1 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.


Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.