Need algorithm for proportional value distributions between weighted buckets

Question

Seems I have problem like Algorithm for constrained proportional distribution but I can't expand solution to my task.

Assume we have X weighted buckets (for example 30000, 17000, 60000) and values Y (for example 2000, 3000, 5000, 8000, 9000, 11000). In matrix all values are marked if they can be distributed to bucket (1 - allowed, 0 - not allowed). We need do distribute values to buckets so the proportions will be as equal as it possible. In this case we will distribute 10000 for 30000 bucket (33%), 8000 for 17000 bucket (47%) and 20000 for 60000 bucket (33%).

My algorithm is:

1. We take values that can be allocated only to 1 bucket and allocate them. Calculate proportion. (5000 for 30000 (16%), 8000 for 17000 (47%) and 11000 for 60000(18%))

2.1. Take values that can be allocated to 2 buckets. Take bucket with min allocation of this two after step 1. Try to fill it with part-values so it will have same proportion as 2nd bucket. Value 3000 can be allocated to 30000 and 17000. At this moment 30000 has 16% and 17000 has 47% so we try to put 3000 to 30000, it will have (5000+3000)/30000 = 26% allocation. Still less than 47%, it's ok. If it will be more than 47% we will fill up to 47% and distribute excess equally between 30000 and 17000 proportionally.

2.2. Take another pair of buckets. And so on. Until we will distribute all 2-bucket values.

3. Take values that can be allocated to 3 buckets. Check min proportion, fill it up to the next. After fill them equally to the 3rd with excess values. Continue until all values will be distributed.

I have feelings that we can get situation where order of distribution (for example 3000->2000 or 2000->3000) can give us different result. So we need to do another steps to avoid that. In this task we don't have any constraints, we only want at the end equal proportions (or as close as possible). Any help with theory for this task? Or may be another algorithm.

Data to check with:
,0,1,2,3,4,5
0,0,9000,9000,4500,4500,1000
1,1000,1,1,1,1,0
2,1000,1,1,1,1,0
3,1000,1,1,1,1,0
4,1000,1,1,1,1,0
5,1000,1,1,1,1,0
6,1000,1,1,1,1,0
7,1000,1,1,1,1,0
8,1000,1,1,1,1,0
9,1000,1,1,1,1,0
10,1000,0,0,0,0,1
and expect to get result:
0.33,0.33,0.17,0.17
,0,1,2,3,4,5
0,0,9000,9000,4500,4500,1000
1,1000,0.33,0.33,0.17,0.17,0
2,1000,0.33,0.33,0.17,0.17,0
3,1000,0.33,0.33,0.17,0.17,0
4,1000,0.33,0.33,0.17,0.17,0
5,1000,0.33,0.33,0.17,0.17,0
6,1000,0.33,0.33,0.17,0.17,0
7,1000,0.33,0.33,0.17,0.17,0
8,1000,0.33,0.33,0.17,0.17,0
9,1000,0.33,0.33,0.17,0.17,0
10,1000,0,0,0,0,1
so we will have one bucket filled for 100%(1000 bucket) and 4 buckets filled equally up to 33%.

Check the algorithm below.

Matrix example
<table>
<tr>
<th>value</th>
<th>30000</th>
<th>17000</th>
<th>60000</th>
</tr>
<tr>
<th>2000</th>
<th>0</th>
<th>1</th>
<th>1</th>
</tr>
<tr>
<th>3000</th>
<th>1</th>
<th>1</th>
<th>0</th>
</tr>
<tr>
<th>5000</th>
<th>1</th>
<th>0</th>
<th>0</th>
</tr>
<tr>
<th>8000</th>
<th>0</th>
<th>1</th>
<th>0</th>
</tr>
<tr>
<th>9000</th>
<th>1</th>
<th>1</th>
<th>1</th>
</tr>
<tr>
<th>11000</th>
<th>0</th>
<th>0</th>
<th>1</th>
</tr>
</table>

Answer 1

The following code is in Python 3 and uses the OR-Tools library.

# reads the CSV export of the .xlsx files, writes CSV to stdout

import csv
import fileinput
import sys

from ortools.linear_solver import pywraplp

# read the input
rows = list(csv.reader(fileinput.input()))
values = [float(row[1]) for row in rows[2:-1]]
buckets = [float(cell) for cell in rows[1][2:-1]]
eligible_member = [[int(cell) for cell in row[2:-1]] for row in rows[2:-1]]

# the solver must have good symmetry breaking, else this formulation is bad
solver = pywraplp.Solver.CreateSolver("solver", "SCIP")

# minimize the max fill ratio minus the min fill ratio
max_fill_ratio = solver.NumVar(0, solver.infinity(), "")
min_fill_ratio = solver.NumVar(0, solver.infinity(), "")
solver.Objective().SetMinimization()
solver.Objective().SetCoefficient(max_fill_ratio, 1)
solver.Objective().SetCoefficient(min_fill_ratio, -1)

# fractional indicator variables -- indicators[i][j] is the fraction of value i
# assigned to bucket j
indicators = [[solver.NumVar(0, cell, "") for cell in row] for row in eligible_member]

# ensure that each value is exactly 100% assigned
for i, value in enumerate(values):
    exact_cover_constraint = solver.Constraint(1, 1)
    for j, bucket in enumerate(buckets):
        exact_cover_constraint.SetCoefficient(indicators[i][j], 1)

# the usual trick for minimizing a maximum -- constrain it to be greater than or
# equal to each value maximized over, let the solver do the rest
for j, bucket in enumerate(buckets):
    max_load_constraint = solver.Constraint(-solver.infinity(), 0)
    max_load_constraint.SetCoefficient(max_fill_ratio, -1)
    min_load_constraint = solver.Constraint(0, solver.infinity())
    min_load_constraint.SetCoefficient(min_fill_ratio, -1)
    for i, value in enumerate(values):
        ratio = value / bucket
        max_load_constraint.SetCoefficient(indicators[i][j], ratio)
        min_load_constraint.SetCoefficient(indicators[i][j], ratio)

# enable to diagnose slowness
if False:
    solver.EnableOutput()

if solver.Solve() != solver.OPTIMAL:
    raise RuntimeError
w = csv.writer(sys.stdout)
w.writerow([solver.Objective().Value()] + buckets)
for i, value in enumerate(values):
    w.writerow([value] + [indicator.solution_value() for indicator in indicators[i]])