How to limit DCs in transportation optimization problem?

Question

In typical transportation optimization problem, we optimize cost for given sources:(Distribution Centers) and demand:(Stores) points.

In below example, how can we add additional constraint to limit no. of Distribution Centers(DC) ?

i.e. Problem :

if dc_lim=10 
total_available_DCs=30
Then Total DCs in optimize network<=dc_lim.

Current Code:

model = ConcreteModel()
model.dual = Suffix(direction=Suffix.IMPORT)

# Step 1: Define index sets
CUS = list(I) # Stores
SRC = list(Supply.keys()) # DCs (Distribution Centres)
K=list(Products)

# Step 2: Define the decision 
model.x = Var(CUS, SRC,K, domain = NonNegativeReals)

# Step 3: Define Objective
model.Cost = Objective(
    expr = sum([T[c,s,k]*model.x[c,s,k] for c in CUS for s in SRC for k in K if (c,s,k) in T]),
    sense = minimize)

# Step 4: Constraints
#Supply Constraint
model.src = ConstraintList()
for s in SRC:
    model.src.add(sum([model.x[c,s,k] for c in CUS for k in K if (c,s,k) in T]) <= Supply[s])

#Demand Constraint
model.dmd = ConstraintList()
for c in CUS:
    for k in K:
        model.dmd.add(sum([model.x[c,s,k] for s in SRC  if (c,s,k) in T]) == Demand[c,k])

        
results = SolverFactory('glpk').solve(model)

Answer 1

It is much better to look at the math instead of a bunch of code. So here we go.

min sum((c,s,k), T(c,s,k)*x(c,s,k))
    sum((c,k), x(c,s,k)) ≤ supply(s)    ∀s
    sum(s, x(c,s,k)) = demand(c,k)      ∀c,k
    x(c,s,k) ≤ M(c,s,k)*y(s)            ∀c,s,k     "if closed don't ship from DC" 
    sum(s,y(s)) ≤ limit                            "limit number of open DCs" 
    y(s) ∈ {0,1}                                   "DC is open or closed"
    x(c,s,k) ≥ 0 
    M(c,s,k) = min(supply(s),demand(c,k))          "constants" 

Transcribing this into code is trivial.