1
$\begingroup$

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.

For example, I have the integer "7" that I need to split according to an array [0.2, 0.3, 0.5] - the result must be [1, 2, 4] - Similarly, if I have the integer "2" that needs to be breakdown according to the array [0.25, 0.25, 0.25, 0.25] the result can be [0, 0, 1, 1] or any other combination until the sum equal 2.

I managed to create the algorithm that is doing exactly the logic described but I was wondering if such a problem cannot be solved with as a linear programming model. Can anyone help me figure out how to produce such model?

I want to express this problem as a LP model because ultimately the breakdown will be more complex and will involve more complex rules. This example is a toy problem.

How can I formulate mathematically this problem so I can encode it into an optimization framework.

Below is the code of the logic :

import numpy as np
from math import isclose, sqrt

def distribute(total, weights):
    scale = float(sum(weights))/total
    return [x / scale for x in weights]

def error_gen(actual, rounded):
    divisor = sqrt(1.0 if actual < 1.0 else actual)
    return abs(rounded - actual) ** 2 / divisor

def round_to_(arr, to):
    if not isclose(sum(arr), to):
        raise ValueError
    n = len(arr)
    rounded = [int(x) for x in arr]
    up_count = to - sum(rounded)
    errors = [(error_gen(arr[i], rounded[i] + 1) - error_gen(arr[i], rounded[i]), i) for i in range(n)]
    rank = sorted(errors)
    for i in range(up_count):
        rounded[rank[i][1]] += 1
    return rounded

if __name__ == "__main__":
    a_case_1 = [.2, .3, .5]
    a_case_2 = [.25, .25, .25, .25]
    a_case_3 = [.25, .25, .5]

    b_case_1 = 7
    b_case_2 = 2
    b_case_3 = 2


    dist_case_1 = distribute(b_case_1, a_case_1)
    print(dist_case_1)
    print(sum(dist_case_1))
    print(round_to_(dist_case_1, b_case_1))

    dist_case_2 = distribute(b_case_2, a_case_2)
    print(dist_case_2)
    print(sum(dist_case_2))
    print(round_to_(dist_case_2, b_case_2))

    dist_case_3 = distribute(b_case_3, a_case_3)
    print(dist_case_3)
    print(sum(dist_case_3))
    print(round_to_(dist_case_3, b_case_3))

------------------ [EDIT] ------------------

Based on the excellent feedback from @LinAlg I did try to implement the suggested formulation with Pyomo however the output doesn't seem to produce something meaningful.

This code is using a different example easier to track. We have a list of cars with some "potentials" and we need to allocate this potential to a specific postal code. the postal code allocation is dictated by some sellout information.

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=NonNegativeIntegers)
model.T = Var([1, 2], domain=NonNegativeReals)


model.objective = Objective(expr=model.T[1] + model.T[2], 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[1]
    )


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[2]
    )


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


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


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


def pyomo_postprocess(options=None, instance=None, results=None):
    model.pprint()
    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)

This code produce the following ouput:

Displaying Solution
--------------------------------------------------------------------------------
5 Set Declarations
    Cars : Dim=0, Dimen=1, Size=3, Domain=None, Ordered=False, Bounds=None
        ['car 1', 'car 2', 'car 3']
    Locations : Dim=0, Dimen=1, Size=3, Domain=None, Ordered=False, Bounds=None
        ['p_code_1', 'p_code_2', 'p_code_3']
    T_index : Dim=0, Dimen=1, Size=2, Domain=None, Ordered=False, Bounds=(1, 2)
        [1, 2]
    X_index : Dim=0, Dimen=2, Size=9, Domain=None, Ordered=False, Bounds=None
        Virtual
    a_index : Dim=0, Dimen=2, Size=9, Domain=None, Ordered=False, Bounds=None
        Virtual

2 Param Declarations
    a : Size=9, Index=a_index, Domain=Any, Default=None, Mutable=False
        Key                   : Value
        ('car 1', 'p_code_1') : 1.4000000000000001
        ('car 1', 'p_code_2') :                2.1
        ('car 1', 'p_code_3') :                3.5
        ('car 2', 'p_code_1') :                0.4
        ('car 2', 'p_code_2') :                0.6
        ('car 2', 'p_code_3') :                1.0
        ('car 3', 'p_code_1') : 2.8000000000000003
        ('car 3', 'p_code_2') :                4.2
        ('car 3', 'p_code_3') :                7.0
    p : Size=3, Index=Cars, Domain=Any, Default=None, Mutable=False
        Key   : Value
        car 1 :     7
        car 2 :     2
        car 3 :    14

2 Var Declarations
    T : Size=2, Index=T_index
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :     0 :   0.0 :  None : False : False : NonNegativeReals
          2 :     0 : 182.0 :  None : False : False : NonNegativeReals
    X : Size=9, Index=X_index
        Key                   : Lower : Value : Upper : Fixed : Stale : Domain
        ('car 1', 'p_code_1') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
        ('car 1', 'p_code_2') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
        ('car 1', 'p_code_3') :     0 :   7.0 :  None : False : False : NonNegativeIntegers
        ('car 2', 'p_code_1') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
        ('car 2', 'p_code_2') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
        ('car 2', 'p_code_3') :     0 :   2.0 :  None : False : False : NonNegativeIntegers
        ('car 3', 'p_code_1') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
        ('car 3', 'p_code_2') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
        ('car 3', 'p_code_3') :     0 :  14.0 :  None : False : False : NonNegativeIntegers

1 Objective Declarations
    objective : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : T[1] + T[2]

3 Constraint Declarations
    sum_maintained : Size=3, Index=Cars, Active=True
        Key   : Lower : Body                                                      : Upper : Active
        car 1 :   7.0 : X[car 1,p_code_1] + X[car 1,p_code_2] + X[car 1,p_code_3] :   7.0 :   True
        car 2 :   2.0 : X[car 2,p_code_1] + X[car 2,p_code_2] + X[car 2,p_code_3] :   2.0 :   True
        car 3 :  14.0 : X[car 3,p_code_1] + X[car 3,p_code_2] + X[car 3,p_code_3] :  14.0 :   True
    t_negative : Size=3, Index=Cars, Active=True
        Key   : Lower : Body                                                                                                    : Upper : Active
        car 1 :  -Inf :  - T[2] - (X[car 1,p_code_1] - 9.8 + X[car 1,p_code_2] - 14.700000000000001 + X[car 1,p_code_3] - 24.5) :   0.0 :   True
        car 2 :  -Inf :                  - T[2] - (X[car 2,p_code_1] - 0.8 + X[car 2,p_code_2] - 1.2 + X[car 2,p_code_3] - 2.0) :   0.0 :   True
        car 3 :  -Inf : - T[2] - (X[car 3,p_code_1] - 39.2 + X[car 3,p_code_2] - 58.800000000000004 + X[car 3,p_code_3] - 98.0) :   0.0 :   True
    t_positive : Size=3, Index=Cars, Active=True
        Key   : Lower : Body                                                                                                : Upper : Active
        car 1 :  -Inf :  X[car 1,p_code_1] - 9.8 + X[car 1,p_code_2] - 14.700000000000001 + X[car 1,p_code_3] - 24.5 - T[1] :   0.0 :   True
        car 2 :  -Inf :                  X[car 2,p_code_1] - 0.8 + X[car 2,p_code_2] - 1.2 + X[car 2,p_code_3] - 2.0 - T[1] :   0.0 :   True
        car 3 :  -Inf : X[car 3,p_code_1] - 39.2 + X[car 3,p_code_2] - 58.800000000000004 + X[car 3,p_code_3] - 98.0 - T[1] :   0.0 :   True

13 Declarations: Cars Locations a_index a p X_index X T_index T objective t_positive t_negative sum_maintained
X : Size=9, Index=X_index
    Key                   : Lower : Value : Upper : Fixed : Stale : Domain
    ('car 1', 'p_code_1') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
    ('car 1', 'p_code_2') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
    ('car 1', 'p_code_3') :     0 :   7.0 :  None : False : False : NonNegativeIntegers
    ('car 2', 'p_code_1') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
    ('car 2', 'p_code_2') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
    ('car 2', 'p_code_3') :     0 :   2.0 :  None : False : False : NonNegativeIntegers
    ('car 3', 'p_code_1') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
    ('car 3', 'p_code_2') :     0 :   0.0 :  None : False : False : NonNegativeIntegers
    ('car 3', 'p_code_3') :     0 :  14.0 :  None : False : False : NonNegativeIntegers
$\endgroup$
1
$\begingroup$

Let $p$ be your integer and let $a \in [0,1]^n$ be your fractions. An example of a formulation is:

\begin{align} \min & \sum_{i=1}^n |x_i-a_i p| \\ \text{s.t.} & \sum_{i=1}^n x_i = p \\ & x \in \mathbb{Z}^n_+ \end{align} The absolute values can be removed with standard methods.

$\endgroup$
  • $\begingroup$ Thanks. May I ask the thought process. I understand how to transform the absolute value to positive and negative constraints however how come did you get the absolute value? $\endgroup$ – Michael Nov 8 '19 at 17:09
  • 1
    $\begingroup$ you want $a_i p$ to be close to $x_i$, and the absolute value is the only distance measure you can easily represent in a MIP (and I just fixed a mistake: $x$ should be a nonnegative integer) $\endgroup$ – LinAlg Nov 8 '19 at 17:34
  • 1
    $\begingroup$ @Michael due to small numerical errors, I do not think you can predict how the branching will go $\endgroup$ – LinAlg Nov 8 '19 at 21:01
  • 1
    $\begingroup$ @Michael seems like you messed up with PositiveIntegers vs NonNegativeIntegers ;) $\endgroup$ – LinAlg Nov 11 '19 at 17:32
  • 1
    $\begingroup$ you need to be more specific in how the fraction should be used as a weight factor $\endgroup$ – LinAlg Nov 11 '19 at 18:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.