Optimization in Python

Solving the candy model

Now let’s see how can we write a program to solve the simple candy model, and we will be using: Pandas, a popular python package for data wrangling, Cplex, an interface to IBM’s cplex solver, and docplex: a modeling package. Solving large scale optimization problem often relies on specialized software, called solvers, and IBM provides a community version of cplex free of charge. If you use pip:

pip install pandas cplex docplex

Or if you use anaconda python distribution:

conda install pandas
conda install -c ibmdecisionoptimization cplex docplex

First, let’s set up the data into Pandas Series and DataFrames.

from docplex.mp.model import Model
import pandas as pd

index = ['Chocolate', 'M&M']
col = ['bundle0', 'bundle1']
A = pd.DataFrame([[2, 4], [3, 2]], index=index, columns=col).T
L = pd.Series([40, 50], index=index)
P = pd.Series([3, 3], index=col)

They look like these:

In [2]: A
Out[2]:
         Chocolate  M&M
bundle0          2    3
bundle1          4    2

In [3]: L
Out[3]:
Chocolate    40
M&M          50
dtype: int64

In [4]: P
Out[4]:
bundle0    3
bundle1    3
dtype: int64

Then we go ahead building the model. The decision we are trying to make, and the unknown we are trying to solve, is for each candy bundle how many we should buy, which should be integer. Therefore we call it decision variable, and in this particular case integer decision variable (line 2). Further more, they has to be non-negative. Therefore we specify some lower bound (lb) for each of the decision variables.
In the end we need to have >=40 chocolate and >=50 M&M. To do that we specify some constraints, using Python generate expression in line 3. Next we need to specify our objective of minimizing cost (line 4).
All done, let’s solve the problem via a .solve call.

There are many guides around showing model construction using Python build-in types such as list, tuple, dictionary etc. Since in real world applications we will use Pandas everyday anyway, we might just as well using Pandas in building model.

mdl = Model('Candies')
X = pd.Series(mdl.integer_var_dict(P.index, lb=0))
C = mdl.add_constraints(mdl.dot(X, A[i]) >= L[i] for i in index)
mdl.minimize(mdl.dot(X, P))
sol = mdl.solve()

Looks like we need to buy 15 bundle0 and 3 bundle1 to reach the cost-minimizing target, while satisfying both of our constraints.

In [5]: sol
Out[5]: docplex.mp.solution.SolveSolution(obj=54,values={_x1:15,_x2:3})

As discussed previously, if the decision variables are allowed to be continuous, we might be able to get even better numerical solution. We can solve the same problem with the integrality removed, which is called solving linear programming (LP) relaxation.

for v in X:
    mdl.set_var_type(v, mdl.continuous_vartype)
    v.set_lb(0)
slp = mdl.solve()

Looks like we now only need to spend 52.5 cents instead of 54 cents from the previous integer solution. However for bundle1 we will buy 2.5 instead of 3.

In [6]: slp
Out[6]: docplex.mp.solution.SolveSolution(obj=52.5,values={_x1:15,_x2:2.5})

_x1 and _x2 in the solution are probably not the most intuitive labels. This one-liner will convert the solution into a nice Pandas Series.

In [7]: pd.Series(sol.get_value_dict(X))
Out[7]:
bundle0    15.0
bundle1     3.0
dtype: float64

Optimization, a geometric intro

Halloween is coming, and I am stockpiling candies. Somehow I have a crystal ball prediction that this year I will need 40 pieces of chocolate and 50 pieces of M&M. In the grocery store there are two kinds of candy bundles, one with 2 chocolate and 3 M&M each, the other one with 4 chocolate and 2 M&M each. They both cost me 3 cents each. Suppose I want to spend as much as possible, and try not to have any left-over candies, how can I figure it out systematically?

This is how we can described it in a mathematical way:

\( Min. 3x + 3y \)
\( s.t. (Subject.to) \)
\( 2x + 4y \geq 40 (chocolate) \)
\( 3x + 2y \geq 50 (M\&M)\)

Well, \(2x + 4y\) geometrically means to project a vector \((x, y)\) to a vector \((2, 4)\). The chocolate constraint basically means that the acceptable solutions must satisfy the condition: once projected to \((2, 4)\), its length needs to be greater or equal to 40. Furthermore, the solutions also need to satisfy a similar constrain for M&M, thus the region of acceptable solutions are in the grey area in the 2nd panel below.

How about the objective? Very similar: project the possible solutions to the objective axis and see which one give us the smallest value.

Well, since we can’t buy a faction of either type of candy bundle, and we have to buy integer numbers of them, the feasible solutions are actually not quite the whole grey area, but the integer grid within it. We call problems with this additional integrality constraint integer problem (IP). Note that the minimal cost we can achieve are slightly different with or without considering integrality.

In a nutshell this is what an optimization problem is all about, and the minimal objective value we can achieve is called optimal solution. In Chinese, solution is 解, which is a angle (角), a cut (刀) and a bull (牛). This is exactly the geometric intuition: rotation some angle draw some axes, make some cuts to determine the feasible region, and finally dissect the bull which is the problem.

⿰ 解
角 ⿱

And is jiě.