code: def flp(I,J,d,M,f,c):
      
  
      model = Model("flp")
  
      x,y = {},{}
      for j in J:
          y[j] = model.addVar(vtype="B", name="y(%s)"%j)
          for i in I:
              x[i,j] = model.addVar(vtype="C", name="x(%s,%s)"%(i,j))
  
      for i in I:
          model.addCons(quicksum(x[i,j] for j in J) == d[i], "Demand(%s)"%i)
  
      for j in M:
          model.addCons(quicksum(x[i,j] for i in I) <= M[j]*y[j], "Capacity(%s)"%i)
  
      for (i,j) in x:
          model.addCons(x[i,j] <= d[i]*y[j], "Strong(%s,%s)"%(i,j))
  
      model.setObjective(
          quicksum(f[j]*y[j] for j in J) +
          quicksum(c[i,j]*x[i,j] for i in I for j in J),
          "minimize")
      model.data = x,y
  
      return model
docs: flp -- model for the capacitated facility location problem
      Parameters:
          - I: set of customers
          - J: set of facilities
          - d[i]: demand for customer i
          - M[j]: capacity of facility j
          - f[j]: fixed cost for using a facility in point j
          - c[i,j]: unit cost of servicing demand point i from facility j
      Returns a model, ready to be solved.