How to use a Pyomo variable result to interpolate a from a list during optimisation?

Question

Let's say you own a fruit store and want to find the optimum number of apples to sell each day. However, there is a different price associated depending on how many apples are sold (you have no say in these bundled price pairs).

e.g.

                                  Day_1                            Day_2
Daily_price_offers = [[(5 apples,$8),(10 apples,$4)] , [(7 apples,$6),(10 apples,$4)], ...]

Now, suppose a customer wishes to purchase 7 apples on Day 1 and 9 apples on Day 2, then they are priced based on the previous bundled price (like a step-wise function). So, in this case they pay $8 on the first day and $6 on the second day.

The Goal

What is the best way to set this up in Pyomo, so that a step-wise function is applied when interpolating? I have created a solution that picks the highest value of the given bundles each day, but I am after the best way for finding the in-between values. Note, this would mean a new variable would need to be defined for the optimal number of apples sold each day.

This example uses integers, but it would be great to be able to handle floating point quantities as well (for different applications).

My Attempt

from pyomo.environ import *

# Inputs
daily_price_offer = [[(5,8),(10,4)],[(7,6),(10,4)],[(4,6),(10,3)]]
max_apples_total = 14
options_per_day = 2

# Define model and solver
model = ConcreteModel()
solver = SolverFactory('gurobi') # Can use gurobi or glpk

# Define variables
model.DAYS = range(len(daily_price_offer))
model.OPTIONS = range(options_per_day)
model.x = Var(model.DAYS, model.OPTIONS, within=Binary)

# Constraint functions
def daily_sale_constraint(model, i):
    "limited to 1 option in a given day"
    return sum(model.x[i,:]) <= 1

# Set constraints
model.weight = Constraint(expr = sum(daily_price_offer[i][j][0] * model.x[i,j] for i in model.DAYS for j in range(options_per_day)) <= max_apples_total)              
model.discharge = Constraint(model.DAYS, rule=daily_sale_constraint)

# Set objective
model.value = Objective(expr = sum(daily_price_offer[i][j][1] * daily_price_offer[i][j][0] * model.x[i,j] for i in model.DAYS for j in range(options_per_day)), sense = maximize)

# Solve
res = solver.solve(model)

# Print results
print('Chosen bundle:',[value(model.x[i,j]) for i in model.DAYS for j in range(options_per_day)])
print('Total revenue: $',model.value())

Answer 1

Interpolation is often described as "piecewise linear functions". There are several ways to do this. The simplest is likely using so-called SOS2 variables (https://en.wikipedia.org/wiki/Special_ordered_set). There are many other ways, however. These approaches can be implemented yourself, or you can use the piecewise library https://pyomo.readthedocs.io/en/stable/library_reference/kernel/piecewise/piecewise.html. See also How to Formulate a Piecewise Step Function in pyomo for an example.

Online help is to be had using:

>>> from pyomo.core import *
>>> help(Piecewise)
Help on class Piecewise in module pyomo.core.base.piecewise:

class Piecewise(pyomo.core.base.block.Block)
 |  Piecewise(*args, **kwds)
 |
 |      Adds piecewise constraints to a Pyomo model for functions of the
 |      form, y = f(x).
 |
 |      Usage:
 |              model.const = Piecewise(index_1,...,index_n,yvar,xvar,**Keywords)
 |              model.const = Piecewise(yvar,xvar,**Keywords)
 |
 |      Keywords:
 |
 |  -pw_pts={},[],()
 |            A dictionary of lists (keys are index set) or a single list
 |            (for the non-indexed case or when an identical set of
 |            breakpoints is used across all indices) defining the set of
 |            domain breakpoints for the piecewise linear
 |            function. **ALWAYS REQUIRED**
 |
 |  -pw_repn=''
 |            Indicates the type of piecewise representation to use. This
 |            can have a major impact on solver performance.
 |            Choices: (Default 'SOS2')
 |
 |               ~ + 'SOS2'      - Standard representation using sos2 constraints
 |               ~   'BIGM_BIN'  - BigM constraints with binary variables.
 |                                 Theoretically tightest M values are automatically
 |                                 determined.
 |               ~   'BIGM_SOS1' - BigM constraints with sos1 variables.
 |                                 Theoretically tightest M values are automatically
 |                                 determined.
 |               ~*+ 'DCC'       - Disaggregated convex combination model
 |               ~*+ 'DLOG'      - Logarithmic disaggregated convex combination model
 |               ~*+ 'CC'        - Convex combination model
 |               ~*+ 'LOG'       - Logarithmic branching convex combination
 |               ~*  'MC'        - Multiple choice model
 |               ~*+ 'INC'       - Incremental (delta) method
 |
 |             + Supports step functions
 |             * Source: "Mixed-Integer Models for Non-separable Piecewise Linear
 |                        Optimization: Unifying framework and Extensions" (Vielma,
 |                        Nemhauser 2008)
 |             ~ Refer to the optional 'force_pw' keyword.
 |
 |  -pw_constr_type=''
 |            Indicates the bound type of the piecewise function.
 |            Choices:
-- More  --

Warning: the BIGM methods are implemented incorrectly. Don't use them.