Pyomo (Transportation theory): including a constraint that compares the value of period 'T' with 'T-1'

I'm trying to solve the classical 'Transportation theory' to optimize the production and transportation of goods from 3 factories to 4 distribution centers, through 2 weeks.

There are:

3 max capacities in each factory (items/week) and each week
3x4 transportation costs from each factory to each distribution center ($/item) and each week
4 weekly demand volume in each distribution center (items/week) and each week

I'm able to solve it with Pyomo:

from pyomo.environ import *
import pyomo.environ as pyo

model = ConcreteModel()

model.T = Set(initialize=[1, 2]) #Weeks

model.i = Set(initialize=['Factory A','Factory B', 'Factory C'])
model.j = Set(initialize=['DC 1', 'DC 2', 'DC 3', 'DC 4'])
model.a = Param(model.T, model.i, initialize={(1, 'Factory A'): 1000, (1, 'Factory B'): 4000, (1, 'Factory C'): 2000, \
    (2, 'Factory A'): 1050, (2, 'Factory B'): 4050, (2, 'Factory C'): 2050}) #Weekly_capacities
model.b = Param(model.T, model.j, initialize={(1, 'DC 1'): 500, (1, 'DC 2'): 900, (1, 'DC 3'): 1800, (1, 'DC 4'):200, \
    (2, 'DC 1'): 600, (2, 'DC 2'): 800, (2, 'DC 3'): 1900, (2, 'DC 4'): 150}) #Weekly_demands

# Transportation costs
costs = {(1, 'Factory A', 'DC 1'): 2,(1, 'Factory A', 'DC 2'): 4, (1, 'Factory A', 'DC 3'): 5, (1, 'Factory A', 'DC 4'): 2, \
    (1, 'Factory B', 'DC 1'): 3, (1, 'Factory B', 'DC 2'): 1, (1, 'Factory B', 'DC 3'): 3, (1, 'Factory B', 'DC 4'): 2, \
    (1, 'Factory C', 'DC 1'): 1, (1, 'Factory C', 'DC 2'): 3, (1, 'Factory C', 'DC 3'): 3, (1, 'Factory C', 'DC 4'): 1,
    (2, 'Factory A', 'DC 1'): 1, (2, 'Factory A', 'DC 2'): 3, (2, 'Factory A', 'DC 3'): 4, (2, 'Factory A', 'DC 4'): 1, \
    (2, 'Factory B', 'DC 1'): 4, (2, 'Factory B', 'DC 2'): 2, (2, 'Factory B', 'DC 3'): 1, (2, 'Factory B', 'DC 4'): 3, \
    (2, 'Factory C', 'DC 1'): 2, (2, 'Factory C', 'DC 2'): 5, (2, 'Factory C', 'DC 3'): 1, (2, 'Factory C', 'DC 4'): 3}

model.d = Param(model.T, model.i, model.j, initialize=costs)

def f_costs(model, T, i, j):
    return model.d[T, i,j]
model.c = Param(model.T, model.i, model.j, initialize=f_costs)

model.x = Var(model.T, model.i, model.j, bounds=(0.0,None)) #Production

def f_supply(model, T, i):
    return sum(model.x[T,i,j] for j in model.j) <= model.a[T,i]
model.supply = Constraint(model.T, model.i, rule=f_supply)

def f_demand(model, T, j):
    return sum(model.x[T,i,j] for i in model.i) >= model.b[T,j]  
model.demand = Constraint(model.T, model.j, rule=f_demand)

def f_objetive(model):
    return sum(model.c[T,i,j]*model.x[T,i,j] for T in model.T for i in model.i for j in model.j)
model.objetive = Objective(rule=f_objetive, sense=minimize)

def pyomo_postprocess(options=None, instance=None, results=None):
    model.x.display()

opt = SolverFactory("glpk")
results = opt.solve(model)
pyomo_postprocess(None, None, results)

Solution:

x : Size=24, Index=x_index
    Key                      : Lower : Value  : Upper : Fixed : Stale : Domain
    (1, 'Factory A', 'DC 1') :   0.0 :    0.0 :  None : False : False :  Reals
    (1, 'Factory A', 'DC 2') :   0.0 :    0.0 :  None : False : False :  Reals
    (1, 'Factory A', 'DC 3') :   0.0 :    0.0 :  None : False : False :  Reals
    (1, 'Factory A', 'DC 4') :   0.0 :    0.0 :  None : False : False :  Reals
    (1, 'Factory B', 'DC 1') :   0.0 :    0.0 :  None : False : False :  Reals
    (1, 'Factory B', 'DC 2') :   0.0 :  900.0 :  None : False : False :  Reals
    (1, 'Factory B', 'DC 3') :   0.0 :  500.0 :  None : False : False :  Reals
    (1, 'Factory B', 'DC 4') :   0.0 :    0.0 :  None : False : False :  Reals
    (1, 'Factory C', 'DC 1') :   0.0 :  500.0 :  None : False : False :  Reals
    (1, 'Factory C', 'DC 2') :   0.0 :    0.0 :  None : False : False :  Reals
    (1, 'Factory C', 'DC 3') :   0.0 : 1300.0 :  None : False : False :  Reals
    (1, 'Factory C', 'DC 4') :   0.0 :  200.0 :  None : False : False :  Reals
    (2, 'Factory A', 'DC 1') :   0.0 :  600.0 :  None : False : False :  Reals
    (2, 'Factory A', 'DC 2') :   0.0 :    0.0 :  None : False : False :  Reals
    (2, 'Factory A', 'DC 3') :   0.0 :    0.0 :  None : False : False :  Reals
    (2, 'Factory A', 'DC 4') :   0.0 :  150.0 :  None : False : False :  Reals
    (2, 'Factory B', 'DC 1') :   0.0 :    0.0 :  None : False : False :  Reals
    (2, 'Factory B', 'DC 2') :   0.0 :  800.0 :  None : False : False :  Reals
    (2, 'Factory B', 'DC 3') :   0.0 : 1900.0 :  None : False : False :  Reals
    (2, 'Factory B', 'DC 4') :   0.0 :    0.0 :  None : False : False :  Reals
    (2, 'Factory C', 'DC 1') :   0.0 :    0.0 :  None : False : False :  Reals
    (2, 'Factory C', 'DC 2') :   0.0 :    0.0 :  None : False : False :  Reals
    (2, 'Factory C', 'DC 3') :   0.0 :    0.0 :  None : False : False :  Reals
    (2, 'Factory C', 'DC 4') :   0.0 :    0.0 :  None : False : False :  Reals

However, I'm struggling to write a constraint that creates an extra cost if the factory must be started up in the second week. This is, in the example 'Factory A' had null production in the first week, but it had 600+150 items in the second week. The model should introduce a fixed cost of $100 for such 'Factory A', thereby modifying the possible outcomes.

I guess that I should create that constraint by comparing the production in period 'T' with period 'T-1', but I'm not able to write it.

Solution

This is all doable. A couple of key ingredients...

You will need 2 new binary variables, indexed appropriately. First, you need to track whether a particular factory is "running" in any period. I did this with a big-M constraint (see code). You will then need a second binary variable to indicate a difference between period t and t-1, which is your "startup" variable.
In pyomo, if the values inside the index set are integers, like time periods, you can "do math" on them, so it would be ok to say model.x[t-1] or something like that. You can also use a couple of neat methods of a pyomo.Set including first() and prev() (there is also last() and next(), but those aren't needed here). This requires the index set is ordered, which is the default and is the order you initialized the set, which, coincidentally, is sorted.
You will need a a constraint that is handles the first period separately, as there is no "look back" from the first period. (see code)
If you are using dictionaries to initialize params (good!), and if the dictionary doesn't cover all combinations, you can include a default= which will cover missing values. (see code, I introduced a 3rd week)

A couple of nits/reminders:

You absolutely must look at the results from the solver and check the status to ensure "optimal" and not something else, or you are looking at junk values
You imported twice, pick your style and stick with it. (I like the pyo. prefix)
Use better names for your sets and pick a convention for capitalization. I like caps for set names, lower for members (see code.)
I usually just number my constraints like C1 or such just because I get confused by the name repetition between variable, function, and constraint name, but to each his own.
I reformatted using black (and then tweaked a bit). I'm still experimenting with that, but I like the way it lines up dictionaries for QA and then I can "fold" them in most IDE/editors

Code:

from pyomo.environ import *
#import pyomo.environ as pyo  # don't import twice...  I prefer this way, but you were using without prefix.

model = ConcreteModel()

# DATA
# Transportation costs
costs = {
    (1, "Factory A", "DC 1"): 2,
    (1, "Factory A", "DC 2"): 4,
    (1, "Factory A", "DC 3"): 5,
    (1, "Factory A", "DC 4"): 2,
    (1, "Factory B", "DC 1"): 3,
    (1, "Factory B", "DC 2"): 1,
    (1, "Factory B", "DC 3"): 3,
    (1, "Factory B", "DC 4"): 2,
    (1, "Factory C", "DC 1"): 1,
    (1, "Factory C", "DC 2"): 3,
    (1, "Factory C", "DC 3"): 3,
    (1, "Factory C", "DC 4"): 1,
    (2, "Factory A", "DC 1"): 1,
    (2, "Factory A", "DC 2"): 3,
    (2, "Factory A", "DC 3"): 4,
    (2, "Factory A", "DC 4"): 1,
    (2, "Factory B", "DC 1"): 4,
    (2, "Factory B", "DC 2"): 2,
    (2, "Factory B", "DC 3"): 1,
    (2, "Factory B", "DC 4"): 3,
    (2, "Factory C", "DC 1"): 2,
    (2, "Factory C", "DC 2"): 5,
    (2, "Factory C", "DC 3"): 1,
    (2, "Factory C", "DC 4"): 3,
    (3, "Factory A", "DC 3"): 2,
    (3, "Factory B", "DC 3"): 1,
}

capacities = {
        (1, "Factory A"): 1000,
        (1, "Factory B"): 4000,
        (1, "Factory C"): 2000,
        (2, "Factory A"): 1050,
        (2, "Factory B"): 4050,
        (2, "Factory C"): 2050,
        (3, "Factory A"): 1000,
        (3, "Factory B"): 500,
    }
demand= {
        (1, "DC 1"): 500,
        (1, "DC 2"): 900,
        (1, "DC 3"): 1800,
        (1, "DC 4"): 200,
        (2, "DC 1"): 600,
        (2, "DC 2"): 800,
        (2, "DC 3"): 1900,
        (2, "DC 4"): 150,
        (3, "DC 3"): 200,
    }

startup_cost = 4
M_prod = max(capacities.values())  # big M value for max production capacity for constraint use

# SETS
model.T = Set(initialize=[1, 2, 3])  # Weeks
model.F = Set(initialize=["Factory A", "Factory B", "Factory C"])  # Factory
model.D = Set(initialize=["DC 1", "DC 2", "DC 3", "DC 4"])         # Destination

# PARAMS
model.cap = Param(model.T, model.F, initialize=capacities, default=0)
model.demand = Param(model.T, model.D, initialize=demand, default=0)
model.cost = Param(model.T, model.F, model.D, initialize=costs, default=10)

# VARS
model.x = Var(model.T, model.F, model.D, bounds=(0.0, None))  # Production
model.running =   Var(model.F, model.T, within=Binary)    # factory f running in period t
model.startup =   Var(model.F, model.T, within=Binary)    # startup cost induced in period t by factory f

# CONSTRAINTS
def f_supply(model, t, f):
    return sum(model.x[t, f, d] for d in model.D) <= model.cap[t, f]
model.supply = Constraint(model.T, model.F, rule=f_supply)


def f_demand(model, t, d):
    return sum(model.x[t, f, d] for f in model.F) >= model.demand[t, d]
model.meet_demand = Constraint(model.T, model.D, rule=f_demand)


# new constraints....
# use a big-M constraint to capture whether a factory is running in period t
def running(model, f, t):
    return sum(model.x[t, f, d] for d in model.D) <= model.running[f, t] * M_prod
model.factory_on = Constraint(model.F, model.T, rule=running)

# link the change in the running variable to the startup cost variable
def startup(model, f, t):
    # handle the first time period seperately....  If it's running, it started in first period
    if t == model.T.first(): 
        return model.startup[f, t] >= model.running[f, t]
    else:
        # look to the previous period's running value to infer startup
        return model.startup[f, t] >= model.running[f, t] - model.running[f, model.T.prev(t)]
model.startup_cost_constraint = Constraint(model.F, model.T, rule=startup)

def f_objetive(model):
    return sum(
        model.cost[t, f, d] * model.x[t, f, d] + startup_cost * model.startup[f, t]
        for t in model.T
        for f in model.F
        for d in model.D
    )
model.objetive = Objective(rule=f_objetive, sense=minimize)


def pyomo_postprocess(options=None, instance=None, results=None):
    print(results)
    for idx in model.x.index_set():
        if model.x[idx]:
            print(f'deliver: {idx}, {value(model.x[idx])}')
    for idx in model.startup.index_set():
        if model.startup[idx]:
            print(f'startup: {idx}')


opt = SolverFactory("glpk")
results = opt.solve(model)
pyomo_postprocess(None, None, results)

Yields:

Problem: 
- Name: unknown
  Lower bound: 11498.0
  Upper bound: 11498.0
  Number of objectives: 1
  Number of constraints: 40
  Number of variables: 55
  Number of nonzeros: 142
  Sense: minimize
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.010598897933959961
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

deliver: (1, 'Factory B', 'DC 2'), 900.0
deliver: (1, 'Factory B', 'DC 3'), 1800.0
deliver: (1, 'Factory C', 'DC 1'), 500.0
deliver: (1, 'Factory C', 'DC 4'), 200.0
deliver: (2, 'Factory A', 'DC 1'), 600.0
deliver: (2, 'Factory A', 'DC 4'), 150.0
deliver: (2, 'Factory B', 'DC 2'), 800.0
deliver: (2, 'Factory B', 'DC 3'), 1900.0
deliver: (3, 'Factory B', 'DC 3'), 200.0
startup: ('Factory A', 2)
startup: ('Factory B', 1)
startup: ('Factory C', 1)
[Finished in 619ms]