Setting a maximum limit on the number of non-zero values in an array of variables

[ConnorGriffin]

I am attempting to create a beer recipe generator using cp-sat that generates recipes based on input parameters and data about different grains. I want the generator to have a set limit on the number of different grains.

For example, I want to generate a recipe with x sweetness using no more than 3 different grain types.

Here is my code right now. It's rough because I'm still learning how to use the library and model my problem, so I apologize.

from ortools.sat.python import cp_model

class SolutionPrinter(cp_model.CpSolverSolutionCallback):
    """Print intermediate solutions."""

    def __init__(self, variables):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__variables = variables
        self.__solution_count = 0

    def OnSolutionCallback(self):
        self.__solution_count += 1
        # for v in self.__variables:
        #     print('%s = %i' % (v, self.Value(v)), end=' ')
        # print()
        for v in self.__variables:
            print('%s = %i' % (v, self.Value(v)), end=' ')
        print()

    def SolutionCount(self):
        return self.__solution_count

target_sweetness_int = 1500

# Grain: Name, Color (SRM * 10000), Potential * 10000, Max Usage Percent, Sweetness
grain_data = [  ('pale 2-row',          18000,      370000,     100,   10),
                ('vienna',              50000,      380000,     100,   15),
                ('munich 10',           100000,     350000,     100,   20),
                ('crystal 40',          400000,     320000,     25,    40),
                ('crystal 120',         1200000,    320000,     15,    100),
                ('chocolate malt (us)', 3000000,    350000,     15,    60)
            ]

num_grains = len(grain_data)
all_grains = range(num_grains)

# Define the model
model = cp_model.CpModel()

# Create the variables for the model
grain_list = [model.NewIntVar(0, 100, 'grain{}'.format(i)) for i in all_grains]

# Grain usage percent must equal 100
model.Add(sum(grain for grain in grain_list) == 100)

# Keep each grain under the max amount
[model.Add(grain_list[grain] <= grain_data[grain][3]) for grain in all_grains]

# Limit sweetness based on target_sweetness
model.Add(
    sum(grain_data[grain][4] * grain_list[grain] for grain in all_grains) == target_sweetness_int
)

# Find all solutions and print the details
solver = cp_model.CpSolver()
solution_printer = SolutionPrinter(grain_list)
status = solver.SearchForAllSolutions(model, solution_printer)
print()
print('Solutions found : %i' % solution_printer.SolutionCount())

This returns 3804 results of grain combinations, but most of them use 4, 5, or 6 different grains.

I don't want to use .Minimize() on the grain count because I need to know all possible solutions, not just one optimal one. This is because I also want to build final beer color into the model, but that requires division against each grain value, and I guess I can't do that with this library. I was planning on feeding all possible solutions into a formula to calculate beer color, and find a best fit that way.

[ghost]

Not sure if it is the best way, but create a "used" flag for each grain and constraint the sum of them:

grain_used = [model.NewBoolVar('grain{}_used'.format(i)) for i in all_grains]
for grain in all_grains:
    model.Add(grain_list[grain] == 0).OnlyEnforceIf(grain_used[grain].Not())
    model.Add(grain_list[grain] > 0).OnlyEnforceIf(grain_used[grain])

model.Add(sum(grain_used) <= 3)

PS: sum(grain for grain in grain_list) reads better as sum(grain_list)

[ConnorGriffin]

Oh wow, that works great. Thanks!

Do you have any ideas on my other problem, regarding calculating the beer color?

Color is calculated like this:

1. Determine the weighted average of the gravity points per pound per gallon (ppg) of the grains based on their potential (grain_data[i][2] and their use percentage (grain_list[i]), and your mash efficiency (leaving this out for now, it's just a percentage).
2. Calculate how many pounds of each grain to add: gravity_points_needed / weighted_ppg * grain_use_percent
3. Calculate each grain's mash_color_units contribution, then add them all together: grain_color * grain_use_pounds / target_volume
4. Apply the Morey Equation to the mash_color_units total: 1.4922 * mash_color_units ^ .6859

I tried building this by calculating variables for each grain like this

# Calculate grain gravity point contributions per pound, factoring in mash efficiency
grain_weighted_sg_vars = [model.NewIntVar(0, 1000000000000, 'grain{}_sg'.format(i)) for i in all_grains]
grain_weighted_sg_list = [model.Add(grain_list[i] * grain_data[i][2] * mash_efficiency_int == grain_weighted_sg_vars[i]) for i in all_grains]

# Calculate how many pounds of each grain to add to hit og - assumes entire beer was that one grain
grain_pounds_to_og_vars = [model.NewIntVar(0, 1000000000000000000, 'grain{}_pounds_to_go'.format(i)) for i in all_grains]
for i in all_grains:
     model.AddDivisionEquality(grain_pounds_to_og_vars[i], sg_needed_int, (grain_list[i] * grain_data[i][2] * mash_efficiency_int))

But this fails with the following error: TypeError: NotSupported: model.GetOrMakeIndex((2701000000 * grain0))

I think I had it written in a list comprehension before instead of a for loop, but it's the same error. Apologies because I am likely missing some variables here, and my scaling is all messed up because I just made the numbers crazy high until I figured out a better scaling method.

[ghost]

create some intermediate variables:

for i in all_grains:
    tmp = model.NewIntVar(0, grain_data[i][3], '')
    model.Add(tmp == (grain_list[i] * grain_data[i][2] * mash_efficiency_int))
    model.AddDivisionEquality(grain_pounds_to_og_vars[i], sg_needed_int, tmp)

[ConnorGriffin]

Ok I will give that a shot. Thank you.

But won't I run into the same problem? model.AddDivisionEquality() won't take a model variable as the third argument. Why does this work but my example does not?

[ghost]

It expects num and denom to be either integer or IntVar, so you have to create intermediate IntVars and set/constraint their values.

[lperron]

Thanks Xian.

Le dim. 29 sept. 2019 à 02:02, Xiang Chen notifications@github.com a écrit :

It expects num and denom to be either integer or IntVar, so you have to create intermediate IntVars and set their values.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/google/or-tools/issues/1603?email_source=notifications&email_token=ACUPL3ILWYMVTAORYPYBUUDQL7V7XA5CNFSM4I3P4QJ2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD73ETOA#issuecomment-536234424, or mute the thread https://github.com/notifications/unsubscribe-auth/ACUPL3LAGA5QDLVJLLBVI2TQL7V7XANCNFSM4I3P4QJQ .

[lperron]

Xiang

Le dim. 29 sept. 2019 à 07:20, Laurent Perron lperron@google.com a écrit :

Thanks Xian.

Le dim. 29 sept. 2019 à 02:02, Xiang Chen notifications@github.com a écrit :

It expects num and denom to be either integer or IntVar, so you have to create intermediate IntVars and set their values.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/google/or-tools/issues/1603?email_source=notifications&email_token=ACUPL3ILWYMVTAORYPYBUUDQL7V7XA5CNFSM4I3P4QJ2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD73ETOA#issuecomment-536234424, or mute the thread https://github.com/notifications/unsubscribe-auth/ACUPL3LAGA5QDLVJLLBVI2TQL7V7XANCNFSM4I3P4QJQ .

[ConnorGriffin]

I feel like a real fool, but I just cannot figure out how to make this work.

I have a slightly reworked/cleaned up version of the code. I am able to calculate the percentage of each grain and the total grain (in pounds) required to hit the desired specific gravity.

from ortools.sat.python import cp_model

# Equipment profile
TARGET_OG = 1.054
MASH_EFFICIENCY = .73
TARGET_VOLUME = 5.75
TARGET_COLOR = 6
TARGET_SWEETNESS = 1
MAX_UNIQUE_GRAINS = 3
SCALING = 100000

class SolutionPrinter(cp_model.CpSolverSolutionCallback):
    """Print intermediate solutions."""

    def __init__(self, variables):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__variables = variables
        self.__solution_count = 0

    def OnSolutionCallback(self):
        self.__solution_count += 1
        # for v in self.__variables:
        #     print('%s = %i' % (v, self.Value(v)), end=' ')
        # print()
        for v in self.__variables:
            print('%s = %i' % (v, self.Value(v)), end=' ')
        print()
        print('Total Weight: {} pounds'.format(self.Value(grain_pounds_total)*100/SCALING))

    def SolutionCount(self):
        return self.__solution_count

# Specify the grains, hardcoded for now
grain_data = [
    {
        "name": "pale 2-row",
        "color": 1.8,
        "potential": 1.037,
        "max_percent": 100,
        "sweetness": 1,
        "bitterness": 0
    },{
        "name": "crystal 40",
        "color": 40,
        "potential": 1.032,
        "max_percent": 25,
        "sweetness": 4,
        "bitterness": 0
    },{
        "name": "crystal 120",
        "color": 120,
        "potential": 1.032,
        "max_percent": 15,
        "sweetness": 10,
        "bitterness": 0
    },{
        "name": "chocolate malt (us)",
        "color": 350,
        "potential": 1.035,
        "max_percent": 15,
        "sweetness": 6,
        "bitterness": 0
    }
]

num_grains = len(grain_data)
all_grains = range(num_grains)
max_ppg = 46 * SCALING

# Scale any non-integer grain properties
grain_data_scaled = []
for grain in grain_data:
    grain['color'] = round(grain['color'] * SCALING)
    grain['potential'] = round((grain['potential'] - 1) * SCALING * MASH_EFFICIENCY)
    #grain['max_percent'] = round(grain['max_percent'] * SCALING)
    grain['sweetness'] = round(grain['sweetness'] * SCALING)
    grain['bitterness'] = round(grain['bitterness'] * SCALING)
    grain_data_scaled.append(grain)

# Scale the equipment profile
equipment = {
    'target_og': int((TARGET_OG - 1) * SCALING),
    'mash_efficiency': int(MASH_EFFICIENCY * SCALING),
    'target_volume': int(TARGET_VOLUME * SCALING),
    'target_sweetness': int(TARGET_SWEETNESS * SCALING),
    'target_color': int(TARGET_COLOR * SCALING)
}
equipment['sg_points_needed'] = int(equipment['target_og'] * equipment['target_volume'])

# Define the model
model = cp_model.CpModel()

# Create the variables for the model
grain_list = [model.NewIntVar(0, 100, 'grain{}'.format(i)) for i in all_grains]

# Grain usage percent must equal 100
model.Add(sum(grain_list) == 100)

# Keep each grain under the max amount
for i in all_grains:
    model.Add(grain_list[i] <= grain_data_scaled[i]['max_percent'])

# Limit the max number of grains to the specified limit
grain_used = [model.NewBoolVar('grain{}_used'.format(i)) for i in all_grains]
for i in all_grains:
    model.Add(grain_list[i] == 0).OnlyEnforceIf(grain_used[i].Not())
    model.Add(grain_list[i] > 0).OnlyEnforceIf(grain_used[i])
model.Add(sum(grain_used) <= MAX_UNIQUE_GRAINS)

# Limit sweetness based on target_sweetness
model.Add(sum(grain_data_scaled[i]['sweetness'] * grain_list[i] for i in all_grains) == equipment['target_sweetness'] * 100)

# Calculate grain gravity point contribution per pound per gallon
grain_weighted_sg_list = [model.NewIntVar(0, max_ppg, 'grain{}_sg'.format(i)) for i in all_grains]
for i in all_grains:
    model.Add(grain_list[i] * grain_data_scaled[i]['potential'] == grain_weighted_sg_list[i])

# Calculate the weighted average ppg
average_ppg = model.NewIntVar(0, max_ppg, 'average_ppg')
model.Add(average_ppg == sum(grain_weighted_sg_list))

# Calculate how many pounds of grain to add total, scaled to 5 decimal places (15.452 = 15452)
# TODO: Find a suitable upper bound
grain_pounds_total = model.NewIntVar(0, 100 * SCALING, 'grain_pounds_total')
model.AddDivisionEquality(grain_pounds_total, equipment['sg_points_needed'], average_ppg)

# Find all solutions and print the details
solver = cp_model.CpSolver()
solution_printer = SolutionPrinter(grain_list)
status = solver.SearchForAllSolutions(model, solution_printer)
print()
print('Solutions found : %i' % solution_printer.SolutionCount())

However now I am getting an error: TypeError: Not an integer linear expression: grain0 when I add the following before the solver. All I'm doing is grain_pounds_total * grain_percent (0-100) for each grain. This is just an intermediary calculation to get to the color calculation. If I can get it working I will move it to a tmp variable as you have suggested.

# Calculate each grain's total weight (will need to divide by 100 later)
grain_pounds_multiplied = [model.NewIntVar(0, 100 * SCALING * 100, 'grain{}_pounds_predivide'.format(i)) for i in all_grains]
for i in all_grains:
    model.Add(grain_pounds_multiplied[i] == grain_pounds_total * grain_list[i])

[ghost]

If you multiply two variables you need to use AddMultiplicationEquality:

for i in all_grains:
    model.AddMultiplicationEquality(grain_pounds_multiplied[i], [grain_pounds_total, grain_list[i]])

Didn't mention it before because I though grain_data[i][2] and mash_efficiency_int where just integers.

[ConnorGriffin]

Thanks again! That did it.

mash_efficency_int was just an int before, but I changed some stuff around which explains why I hadn't run into this before. I had been multiplying inside of model.Add() all day, but only ever using a single IntVar at a time.

I really appreciate the help. May I keep this issue open for a little while longer should I run into other roadblocks in the next couple days?