import pyomo.environ as pyo
model = pyo.ConcreteModel(name='Slick Oil')

# Decision variables
model.w1_rA = pyo.Var(within=pyo.NonNegativeReals)
model.w1_rB = pyo.Var(within=pyo.NonNegativeReals)

model.w2_rB = pyo.Var(within=pyo.NonNegativeReals)
model.w2_rC = pyo.Var(within=pyo.NonNegativeReals)
model.w2_rD = pyo.Var(within=pyo.NonNegativeReals)

model.w3_rB = pyo.Var(within=pyo.NonNegativeReals)
model.w3_rC = pyo.Var(within=pyo.NonNegativeReals)
model.w3_rE = pyo.Var(within=pyo.NonNegativeReals)

model.w4_rA = pyo.Var(within=pyo.NonNegativeReals)
model.w4_rE = pyo.Var(within=pyo.NonNegativeReals)

model.w5_rD = pyo.Var(within=pyo.NonNegativeReals)
model.w5_rE = pyo.Var(within=pyo.NonNegativeReals)

model.w6_rC = pyo.Var(within=pyo.NonNegativeReals)
model.w6_rE = pyo.Var(within=pyo.NonNegativeReals)

# Objective function
model.obj = pyo.Objective(expr=7*model.w1_rA + 5*model.w1_rB
                          + 14*model.w2_rB + 18*model.w2_rC + 13*model.w2_rD
                          + 6*model.w3_rB + 10*model.w3_rC + 12*model.w3_rE
                          + 10*model.w4_rA + 14*model.w4_rE
                          + 11*model.w5_rD + 18*model.w5_rE
                          + 8.5*model.w6_rC + 10.5*model.w6_rE,
                          sense=pyo.minimize)

# Well constraints
model.w1 = pyo.Constraint(expr=model.w1_rA + model.w1_rB <= 40)
model.w2 = pyo.Constraint(expr=model.w2_rB + model.w2_rC + model.w2_rD <= 50)
model.w3 = pyo.Constraint(expr=model.w3_rB + model.w3_rC + model.w3_rE <= 40)
model.w4 = pyo.Constraint(expr=model.w4_rA + model.w4_rE <= 100)
model.w5 = pyo.Constraint(expr=model.w5_rD + model.w5_rE <= 100)
model.w6 = pyo.Constraint(expr=model.w6_rC + model.w6_rE <= 100)

# Refinery constraints
model.r1 = pyo.Constraint(expr=model.w1_rA + model.w4_rA <= 100)
model.r2 = pyo.Constraint(expr=model.w1_rB + model.w2_rB + model.w3_rB <= 40)
model.r3 = pyo.Constraint(expr=model.w2_rC + model.w3_rC + model.w6_rC <= 80)
model.r4 = pyo.Constraint(expr=model.w2_rD + model.w5_rD <= 100)
model.r5 = pyo.Constraint(expr=model.w3_rE + model.w4_rE + model.w5_rE + model.w6_rE <= 80)

# Demand constraint
model.demand = pyo.Constraint(expr=model.w1_rA + model.w1_rB
                              + model.w2_rB + model.w2_rC + model.w2_rD
                              + model.w3_rB + model.w3_rC + model.w3_rE
                              + model.w4_rA + model.w4_rE
                              + model.w5_rD + model.w5_rE
                              + model.w6_rC + model.w6_rE == 100)

model.write("slick_oil.mps", io_options={"symbolic_solver_labels": True})