Convert absolute value in pyomo

[Michaël Lambé]

I am working on a project where I need to breakdown an integer value according to an array of percentage values. My end array must contain integer value and the sum of the array must be equal to the initial integer.

We assume p be our integer and let 𝑎 be our fractions. An example of a formulation is:

mins.t.∑𝑖=1𝑛|𝑥𝑖−𝑎𝑖𝑝|∑𝑖=1𝑛𝑥𝑖=𝑝𝑥∈ℤ𝑛+

I am trying to build such logic using pyomo. Below is my code but I cannot get my head around the conversion of the absolute value.

In this example, I have a series of cars for which I have some potentials sales (integer value) and I need to distribute these sales to postal code level according to some sellouts.

The function `distribute` make sure the sellout is scaled according to the potentials.

I would really appreciate your support. I am relatively new to PyOmo and very keen to learn. Thanks in advance.

from pyomo.environ import *


def distribute(total, weights):
    scale = float(sum(weights.values())) / total
    return {k: v / scale for k, v in weights.items()}


Cars = ["car 1", "car 2", "car 3"]
Locations = ["p_code_1", "p_code_2", "p_code_3"]
POTENTIALS = {"car 1": 7, "car 2": 2, "car 3": 14}
SELLOUTS = {"p_code_1": 0.2, "p_code_2": 0.3, "p_code_3": 0.5}

SELLOUTS_PER_P_CODE = {}

for car in Cars:
    pot = POTENTIALS[car]
    scaled_sellout = distribute(pot, SELLOUTS)
    t = {(car, p_code): v for p_code, v in scaled_sellout.items()}
    SELLOUTS_PER_P_CODE.update(t)

model = ConcreteModel(name="Breakdown Potential to Post Code")

model.Cars = Set(initialize=Cars)
model.Locations = Set(initialize=Locations)

model.a = Param(model.Cars, model.Locations, initialize=SELLOUTS_PER_P_CODE)
model.p = Param(model.Cars, initialize=POTENTIALS)

model.X = Var(model.Cars, model.Locations, within=PositiveIntegers)
model.T = Var(model.Cars, within=PositiveIntegers)


def objective_rule(model):
    return sum(model.T[i] for i in model.Cars)


model.objective = Objective(rule=objective_rule, sense=minimize)


def t_positive_rule(model, i):
    return (
        sum(model.X[i, j] - model.a[i, j] * model.p[i] for j in model.Locations)
        <= model.T[i]
    )


model.t_positive = Constraint(model.Cars, rule=t_positive_rule)


def t_negative_rule(model, i):
    return (
        sum(model.X[i, j] - model.a[i, j] * model.p[i] for j in model.Locations)
        >= -model.T[i]
    )


model.t_negative_rule = Constraint(model.Cars, rule=t_negative_rule)


def sum_maintained_rule(model, j):
    return sum(model.X[j, i] for i in model.Locations) == model.p[j]


model.sum_maintained = Constraint(model.Cars, rule=sum_maintained_rule)


def pyomo_postprocess(options=None, instance=None, results=None):
    print("Total Objective:", model.objective())
    model.X.display()


if __name__ == "__main__":
    opt = SolverFactory("glpk")
    results = opt.solve(model)
    results.write()
    print("\nDisplaying Solution\n" + "-" * 80)
    pyomo_postprocess(None, model, results)

[Michaël Lambé]

The formulation doesn't appear to be good, here is a picture:

[Vincent Reinbold]

HI Michaël,

I am not sure it is what you want, but minimizing an absolute value is usually a MILP.

You need to consider two new positive variables (the positive and the negative part of what is inside the absolute value), and a binary variable such as u=1 if positive and u=0 otherwise.

This linearization trick can be done using the big M method, such as :

min |y|    =>   min (p + n), st.  y = p - n, p <= u.M, n<=(1-u).M , p, n \in R+ , u \in {0;1} and M a big parameter.

try to adapt this to you problem (with indexes and all), hope it can help,

Greetings

[Erling D. Andersen]

Why not do

min sum_i t_i

st.   x_i - a_i*p <= t_i

       -(x_i - a_i*p) <= t_i 

       other constraints