MIP with piecewise functions

MIP with piecewise functions.

Sometimes a relationship in optimization is non-linear and makes it unsolvable especially for non-commercial solver.

A common trick is converting the non-linear relationship to a piecewise function, so that for each interval of the function it is a linear relationship. It would need to model for each interval separately with extra constraints. The more intervals used for the piece-wise function, the better the approximation is, but obviously more constraints mean.

An example: A farmer has a block of land and chooses how much land to allocate to growing corn and wheat to maximize profit, subject to land and budget constraints. The cost per hectare decreases as more is planted (economies of scale) so the cost function is not linear and can be divided into a piecewise function.

Problem Setup

Decision variables:

•         x: hectares of corn

•         y: hectares of wheat

Constraints:

•         Total land available: x+y <= 40

•         Total cost of planting corn and wheat must not exceed budget $35,000

Piecewise cost functions:

•         Corn:

As x increases, the cost increase slows down. The portion of x within the below buckets increases cost differently.

            x in [0,10]: cost increase = 1000x

            x in (10,20]: cost increase = 800x

            x in (20,30]: cost increase = 600x

E.g. if x = 15, the cost = 1000 * 10 + 800* 5, i.e. cost = 800 * x + 2000

  if x = 25, the cost = 1000 * 10 + 800* 10 + 5 * 600, i.e. cost = 600*x + 6000

•         Wheat:

            y in [0,10]: cost = 900y

            y in (10,20]: cost = 700y

            y in (20,30]: cost = 500y

•         profits:

            Corn: profit = $1500/hectare

            Wheat: profit = $1700/hectare

To solve the problem, introduce additional variables:

            zx1, zx2, zx3 : binary indicators for indicating which interval x falls in.

            zy1, zy2, zy3 : binary indicators for indicating which interval y falls in.

Note x or y can fall in only one interval. If zx2 is 1, then zx1 and zx3 must be 0.

If x falls in zx2, it also needs to enforce constraints x in (10,20]. This is done using the big M trick.

Constraints:

•         x + y <= 40 # Land constraint

Cost for x

•         zx1 + zx2 + zx3 = 1 # fall in only one interval

•         x >= zx1 * 0

•         x >= zx2 * 10

•         x >= zx3 * 20

•         x <= zx1 * 10 + M * (1-zx1)  # if zx1 is chosen, x must be between 0 and 10.

•         x <= zx2 * 20 + M * (1-zx2)

•         x <= zx3 * 30 + M * (1-zx3)

•         cost_x = 1000 * x * zx1 + (800 * x + 2000) * zx2 + (600 * x + 6000) * zx3 #quadratic

Cost for y

•         zy1 + zy2 + zy3 = 1 # fall in only one interval

•         y >= zy1 * 0

•         y >= zy2 * 10

•         y >= zy3 * 20

•         y <= zy1 * 10 + M * (1-zy1)  # if zy1 is chosen, y must be between 0 and 10.

•         y <= zy2 * 20 + M * (1-zy2)

•         y <= zy3 * 30 + M * (1-zy3)

•         cost_y = 900 * y * zy1 + (700 * y + 1500) * zy2 + (500 * y + 4500) * zy3 # quadratic

Cost constraint

•         cost_x + cost_y <= $35,000

Objective: Maximize profit

•         Profit=1500x + 1700y

!!NOTE the above solution comes with quadratic constraints for the costs of x and y.

The quadratic constraints cannot be handled by many solvers like in ORTOOLs. However it can be easily handled by Gurobi.

Here is a script in Gurobi. A linear solution is covered later.

from gurobipy import Model, GRB

# Parameters

L = 40  # Total land available

B = 35000   # Budget

M = 1000 # Big M

# Create model

model = Model("FarmOptimization")

# Decision variables

x = model.addVar(lb=0, ub=L, name="x")  # hectares of corn

y = model.addVar(lb=0, ub=L, name="y")  # hectares of wheat

# Binary variables for piecewise intervals

zx1 = model.addVar(vtype=GRB.BINARY, name="zx1")

zx2 = model.addVar(vtype=GRB.BINARY, name="zx2")

zx3 = model.addVar(vtype=GRB.BINARY, name="zx3")

zy1 = model.addVar(vtype=GRB.BINARY, name="zy1")

zy2 = model.addVar(vtype=GRB.BINARY, name="zy2")

zy3 = model.addVar(vtype=GRB.BINARY, name="zy3")

# Cost variables

cost_x = model.addVar(lb=0, name="cost_x")

cost_y = model.addVar(lb=0, name="cost_y")

# Land constraint

model.addConstr(x + y <= L, "LandConstraint")

# Interval selection constraints

model.addConstr(zx1 + zx2 + zx3 == 1, "CornInterval")

model.addConstr(zy1 + zy2 + zy3 == 1, "WheatInterval")

# Interval bounds for x (corn)

model.addConstr(x >= 0 * zx1)

model.addConstr(x >= 10 * zx2)

model.addConstr(x >= 20 * zx3)

model.addConstr(x <= 10 + M * (1 - zx1))

model.addConstr(x <= 20 + M * (1 - zx2))

model.addConstr(x <= 30 + M * (1 - zx3))

# Interval bounds for y (wheat)

model.addConstr(y >= 0 * zy1)

model.addConstr(y >= 10 * zy2)

model.addConstr(y >= 20 * zy3)

model.addConstr(y <= 10 + M * (1 - zy1))

model.addConstr(y <= 20 + M * (1 - zy2))

model.addConstr(y <= 30 + M * (1 - zy3))

# Cost expressions (quadratic terms allowed in Gurobi)

model.addConstr(

cost_x == 1000 * x * zx1 + (800 * x + 2000) * zx2 + (600 * x + 6000) * zx3,

"CostCorn"

)

model.addConstr(

cost_y == 900 * y * zy1 + (700 * y + 1500) * zy2 + (500 * y + 4500) * zy3,

"CostWheat"

)

# Budget constraint

model.addConstr(cost_x + cost_y <= B, "Budget")

# Objective: Maximize profit

model.setObjective(1500 * x + 1700 * y, GRB.MAXIMIZE)

# Solve

model.optimize()

# Output results

if model.status == GRB.OPTIMAL:

print(f"Optimal solution found:")

print(f"Corn hectares (x): {x.X}")

print(f"Wheat hectares (y): {y.X}")

print(f"Profit: {model.ObjVal}")

print(f"Cost for corn: {cost_x.X}")

print(f"Cost for wheat: {cost_y.X}")

else:

print("No optimal solution found.")

Optimal solution found (tolerance 1.00e-04)

Best objective 6.600000000000e+04, best bound 6.600000000000e+04, gap 0.0000%

Optimal solution found:

Corn hectares (x): 10.0

Wheat hectares (y): 30.0

Profit: 66000.0

Cost for corn: 10000.0

Cost for wheat: 19500.0

Linearize the solution

To replace the quadratic constraints, it needs to rewrite the variables and constraints.

Introduce additional variables:

            x1, x2, x3: the x value portion for each interval

            zx1, zx2, zx3 : binary indicator for using the interval

            y1, y2, y3: the y value portion for each interval

            zy1, zy2, zy3 : binary indicator for using the interval

Note x1 + x2 + x3 = x, and x1 must be fully used before using x2, and so on. Therefore, zx1 must be 1 if zx2 is 1.

Constraints:

•         x + y <= 35  #land constraint

•         x1 + x2 + x3 = x

•         y1 + y2 + y3 = y

Cost for x

•         0<= x1 <= 10

•         0<= x2 <= 10

•         0<= x3 <= 10

•         zx3 <= zx2  # if an interval is used, earlier intervals must be used as well

•         zx2 <= zx1

•         x1 >= 10 * zx2  # if an interval is used, earlier intervals must be used FULLY

•         x2 >= 10 * zx3

•         cost_x = 1000 * x1 + 800 * x2 + 600 * x3

Cost for y

•         0<= y1 <= 10

•         0<= y2 <= 10

•         0<= y3 <= 10

•         zy3 <= zy2

•         zy2 <= zy1

•         y1 >= 10 * zy2

•         y2 >= 10 * zy3

•         cost_y = 900 * y1 + 700 * y2 + 500 * y3

Cost constraint

•         cost_x + cost_y <= 35,000

Objective: Maximize profit

•         Profit = 1500x + 1700y

Now this can be solved by ORTools because its all linear.

from ortools.linear_solver import pywraplp

# Create solver

solver = pywraplp.Solver.CreateSolver('SCIP')

# Parameters

L = 40

B = 35000

# Decision variables

x = solver.NumVar(0, L, 'x')

y = solver.NumVar(0, L, 'y')

# Interval variables

x1 = solver.NumVar(0, 10, 'x1')

x2 = solver.NumVar(0, 10, 'x2')

x3 = solver.NumVar(0, 10, 'x3')

y1 = solver.NumVar(0, 10, 'y1')

y2 = solver.NumVar(0, 10, 'y2')

y3 = solver.NumVar(0, 10, 'y3')

# Binary indicators

zx1 = solver.BoolVar('zx1')

zx2 = solver.BoolVar('zx2')

zx3 = solver.BoolVar('zx3')

zy1 = solver.BoolVar('zy1')

zy2 = solver.BoolVar('zy2')

zy3 = solver.BoolVar('zy3')

# Total land constraint

solver.Add(x + y <= L)

# Interval sum constraints

solver.Add(x1 + x2 + x3 == x)

solver.Add(y1 + y2 + y3 == y)

# Interval bounds and sequencing for corn

solver.Add(x1 <= 10)

solver.Add(x2 <= 10)

solver.Add(x3 <= 10)

solver.Add(zx3 <= zx2)

solver.Add(zx2 <= zx1)

solver.Add(x1 >= 10 * zx2)

solver.Add(x2 >= 10 * zx3)

# Interval bounds and sequencing for wheat

solver.Add(y1 <= 10)

solver.Add(y2 <= 10)

solver.Add(y3 <= 10)

solver.Add(zy3 <= zy2)

solver.Add(zy2 <= zy1)

solver.Add(y1 >= 10 * zy2)

solver.Add(y2 >= 10 * zy3)

# Cost expressions

cost_x = solver.NumVar(0, solver.infinity(), 'cost_x')

cost_y = solver.NumVar(0, solver.infinity(), 'cost_y')

solver.Add(cost_x == 1000 * x1 + 800 * x2 + 600 * x3)

solver.Add(cost_y == 900 * y1 + 700 * y2 + 500 * y3)

# Budget constraint

solver.Add(cost_x + cost_y <= B)

# Objective: Maximize profit

profit = 1500 * x + 1700 * y

solver.Maximize(profit)

# Solve

status = solver.Solve()

# Output

if status == pywraplp.Solver.OPTIMAL:

print("Optimal solution found:")

print(f"x (corn): {x.solution_value()}")

print(f"y (wheat): {y.solution_value()}")

print(f"x1: {x1.solution_value()}, x2: {x2.solution_value()}, x3: {x3.solution_value()}")

print(f"y1: {y1.solution_value()}, y2: {y2.solution_value()}, y3: {y3.solution_value()}")

print(f"Cost_x: {cost_x.solution_value()}, Cost_y: {cost_y.solution_value()}")

print(f"Profit: {solver.Objective().Value()}")

else:

print("No optimal solution found.")