An answer to this question on Stack Overflow.

Question

Consider this example problem:

Suppose the cost for setting up a factory to generate a pencil is 1000 and to generate a pen is 2000. The profit for each pencil is 10 and the profit for each pen is 12. For each month, the factory can generate 100 items (pen or pencil). How can I include the initial cost in order to maximize the profit of the factory?

I can think of:

Maximize 10x1 + 12x2 - (1000b1 + 2000b2)
x1+x2 <= 100

Where b1 and b2 are binary variables indicating whether to generate a pen or pencil. But I don't know how to relate them to x1 and x2, so that, for example, if x1 >= 0 then b1 =1, otherwise b1=0

Answer

This is a fixed-cost integer program.

Let's introduce the following variables:

• x: The continuous number of units of the good produced
• y: A zero-one variable indicating whether or not fixed costs are incurred
• c: The per unit revenue of producing X
• f: The fixed cost incurred from producing a nonzero quantity of X regardless of how many units are produced
• m: A large number

The cost function with fixed costs is is:

[![Fixed cost cost function][1]][1]

We transform this into the following program:

max cx - fy
s.t.   x<=my
       0<=x
       y={0,1}

Note that if x=0, then the cost function will set y=0 in order to avoid paying f. However, if x>0, then the x<=my constraint forces y=1. The only issue is that m must be an upper bound on x. We can achieve this either through carefully thinking about what values x can take, or setting m to be a very large number.

We can solve the program using cvxpy as follows. Note that if x1+x2<=100, then it is better to produce nothing. Since this is boring, I've increased the limit from 100 to 300.

#!/usr/bin/env python3
import cvxpy as cp
b1 = cp.Variable(boolean=True)
b2 = cp.Variable(boolean=True)
x1 = cp.Variable()
x2 = cp.Variable()
m = 100000 #A very large number
obj = 10*x1+12*x2 - 1000*b1 - 2000*b2
cons = []
cons += [x1>=0]
cons += [x2>=0]
cons += [x1+x2<=300]
cons += [x1<=m*b1]
cons += [x2<=m*b2]
prob = cp.Problem(cp.Maximize(obj), cons)
objval = prob.solve()
print("Objective value: {0}".format(objval))
print("x1 = {0}".format(x1.value))
print("x2 = {0}".format(x2.value))
print("b1 = {0}".format(b1.value))
print("b2 = {0}".format(b2.value))

[1]: https://i.sstatic.net/bMYY2.png