from pyomo.environ import *
from pyomo.opt import *
opt = solvers.SolverFactory("glpk") 

Jobs  = [0, 1, 2, 3, 4]
Slots = [0, 1, 2, 3, 4]

Due     = [2, 2, 3, 1, 3]
Penalty = [19, 12, 34, 30, 22]

model = ConcreteModel()

model.a = Var(Jobs, Slots, within=Boolean)

model.z = Objective(
    expr = sum(model.a[j,s]*(s-Due[j])*Penalty[j]
               for j in Jobs for s in Slots if s>Due[j]), 
    sense=minimize)

def all_jobs_assigned_rule (model, j):
    return sum(model.a[j,s] for s in Slots) == 1

model.jobs = Constraint(Jobs, rule=all_jobs_assigned_rule)

def all_slots_assigned_rule (model, s):
    return sum(model.a[j,s] for j in Jobs) == 1

model.slots = Constraint(Slots, rule=all_slots_assigned_rule)

results = opt.solve(model)

model.a.get_values()

{(0, 0): 1.0,
 (0, 1): 0.0,
 (0, 2): 0.0,
 (0, 3): 0.0,
 (0, 4): 0.0,
 (1, 0): 0.0,
 (1, 1): 0.0,
 (1, 2): 1.0,
 (1, 3): 0.0,
 (1, 4): 0.0,
 (2, 0): 0.0,
 (2, 1): 0.0,
 (2, 2): 0.0,
 (2, 3): 1.0,
 (2, 4): 0.0,
 (3, 0): 0.0,
 (3, 1): 1.0,
 (3, 2): 0.0,
 (3, 3): 0.0,
 (3, 4): 0.0,
 (4, 0): 0.0,
 (4, 1): 0.0,
 (4, 2): 0.0,
 (4, 3): 0.0,
 (4, 4): 1.0}

model.z.expr()

22.0