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

P = ['P1', 'P2', 'P3']
M = ['M1', 'M2', 'M3']

s = {'P1':3, 'P2':5, 'P3':1}
d = {'M1':2, 'M2':2, 'M3':5}

c = {('P1','M1'):10, ('P1','M2'):5,  ('P1','M3'):15,
     ('P2','M1'):15, ('P2','M2'):10, ('P2','M3'):20,
     ('P3','M1'):5,  ('P3','M2'):5,  ('P3','M3'):15}

model = ConcreteModel()
model.x = Var(P, M, within=NonNegativeReals)

model.z = Objective(expr = sum(c[i,j]*model.x[i,j] for i in P for j in M), sense=minimize)

def supply_rule (model, i):
    return sum(model.x[i,j] for j in M) <= s[i]
model.supply = Constraint(P, rule=supply_rule)

def demand_rule (model, j):
    return sum(model.x[i,j] for i in P) >= d[j]
model.demand = Constraint(M, rule=demand_rule)

model.dual = Suffix(direction=Suffix.IMPORT)
results = opt.solve(model)
model.x.get_values()

{('P1', 'M1'): 0.4634679251867636,
 ('P1', 'M2'): 0.8561172866483512,
 ('P1', 'M3'): 1.6804148178139067,
 ('P2', 'M1'): 0.5365320260507703,
 ('P2', 'M2'): 1.1438827029955108,
 ('P2', 'M3'): 3.3195851417491173,
 ('P3', 'M1'): 1.0000000288526543,
 ('P3', 'M2'): 0.0,
 ('P3', 'M3'): 0.0}

model.z.expr()

124.99999835228026

for j in M:
    print(j, model.dual[model.demand[j]])

M1 15.043624766149037
M2 10.043624768629797
M3 20.043624770065662

for i in P:
    print(i, model.dual[model.supply[i]])

P1 -5.043624771556779
P2 -0.043624770820542835
P3 -10.043624768654922