Linear Programming for Data Science

Original article was published by Yudhiesh on Deep Learning on Medium


Linear Programming for Data Science

As someone aspiring to be a Data Scientist I overlooked one of the most basic and important topic that is needed to become a Data Scientist which is Linear Programming. Linear Programming is simply a technique to depict complex relationships through linear functions and then find the optimum points or in layman’s terms performing optimisation of a problem. Optimisation is part of our daily lives as we have to live with a finite amount of resources as well as time and we need to make the most of it in order to survive. Optimisation could come in simple forms such as planning our shopping trip beforehand all the way up to creating an algorithm to connect airline passengers to their destination with all the travel restrictions during the current COVID-19 pandemic. Here I will go over a problem to help you better understand how to use Linear Programming for any problem.

NOTE: I am assuming that you are familiar with high school mathematics as I will not be going into detail explaining certain parts.

Github repository with more examples

The Snowball Problem

How big would a snowball be initially in order to knock down a tree after rolling down a hill after 30 seconds? Let’s assume it takes 25000N to knock down a tree.

Formula’s and parameters:

Here are the formula’s and the parameters we will use to help us solve this problem.

Parameters:

Initial snowball conditions:

Formula for radius:

Function for dynamic snowball:

Function for objective function:

I know it looks quite intimidating at first being thrown a bunch of formulas and parameters but bear with me on this, it will become easier to understand once we start solving it.

Let’s start coding!

Here are the imports needed for this. NumPy is a python library used for scientific computing. SciPy is a python library used for mathematical functions.

import numpy as np
import
pulp as p
from scipy.optimize import minimize
from
scipy.integrate import odeint

Decision Variables and Constraints

Here we are defining the parameters needed or the decision variables. These variables will decide the output. The constraints are the restrictions or limitations on the decision variables. They act as a limiter to limit the value of the decision variables. Here the constrains are the Target Force and the Time.

# define the parameters
K0 = 85
beta = 0.07
C_d = 0.3
g = 9.8
rho = 350
theta = np.radians(5)
rho_a = 0.9

# initial snowball conditions
m0 = 10
v0 = 0
r0 = (m0/(4/3.0* np.pi *rho))**(1/3.0)
s0 = 0

# Target Force
F_d = 25000

# Set up time array to solve for 30 seconds
t = np.linspace(0,30)

Here we are defining the function of the equation of motion of a snowball rolling down the hill. We unpack the variables mass, radius, distance moved by the object as well as the velocity from w.

We also unpack the parameters from p and use them in f which stands for the functions of the four differential equations. The functions are stored in a list which we can index in order to access later on. While it looks complicated all those equations are just the code for the differential equations and you have to be extra careful when inputting them as it is quite a tedious task and easy to make a mistake.

# define the function of dynamics snowball, the equations of motion
def snowball_dynamics(w,t,p):

#Unpack state variables
M, r, s, v = w

#Unpack parameters
K0, C_d, g, rho, theta, rho_a, beta = p

#Make an array of the right hand sides of the four differential equations that make up our system
f = [beta * K0 * np.exp(-beta*t),
(beta * K0 * np.exp(-beta*t))/(4*np.pi*rho*r**2),
v,
(-15*rho_a*C_d)/(56*rho)*1/r*v**2-23/7*1/r*beta*K0*np.exp(-beta*t)/(4*np.pi*rho*r**2)*v+5/7*g*np.sin(theta)]
return f

The objective function and Non-negative restriction

The objective function is used to set the objective of making decisions. In this example, we are trying to figure out or minimise the size of a snowball that would knock down a tree. Big is a bit subjective so we instead will figure out the mass and radius of the snowball. When we talk about minimising a function what I mean is what is the minimum possible value of that function to satisfy the question? When dealing with Gradient Descent the objective is to minimise the cost function in order to get the lowest possible error of the model on the training instances.

Non-negative restrictions ensure that the decision variables should always take non-negative values or that they should be greater than or equal to 0. This is done in the second last line where obj is declared.

#Define the objective function, minimise the output of this function by changing the initial snowball mass
def objective(m0):
#load parameters
p = [m0, C_d, g, rho, theta, rho_a, beta]
#get the initial radius from initial mass
r0 = (m0/(4/3.0*np.pi*rho))**(1/3.0)
#set initial guesses
w0 = [m0, r0, s0, v0]
#integrate forward for 30 seconds
sol = odeint(snowball_dynamics, w0, t, args=(p,))
#calculate kinetic energy at the end of the run
#
sol[:,0][-1]
# we are accessing the last values from the first column which stores mass
# column 3 contains the values of velocity
ke = 0.5 * sol[:,0][-1]* sol[:,3][-1]**2 #calculate force required to stop snowball in one snowball radius
F = ke/ sol[:,1][-1]
#Compare to desired force: This should be equal to zero when we are done
obj = (F - F_d)**2
return obj

Finally run the function to get the results

Here we are minimising the objective function with the parameter m0. As we are looking for the initial mass we take the first value in the result and use it to calculate the initial mass in kilograms and pounds as well as the initial radius in the form of centimetres and inches.

# Call optimisation using the function defined above
res = minimize(objective, m0, options ={"disp": True})

# Get the optimized initial mass from solution
m0_opt = res.x[0]

# calculate the optimised initial radius from the initial mass
r0_opt = (m0_opt/(4/3.0*np.pi*rho))**(1/3.0)

print(f'Initial mass : {m0_opt}kg {m0_opt*2.02}lbs')
print(f'Initial radius : {r0_opt*100}cm {r0_opt*39.37}inches')

Output:

There is our answer of an initial mass of roughly 54.5kgs and an initial radius of 33.4cm.

Optimization terminated successfully.
Current function value: 0.000000
Iterations: 9
Function evaluations: 33
Gradient evaluations: 11
Initial mass : 54.4727546940319kg 110.03496448194444lbs
Initial radius : 33.36882350231841cm 13.137305812862758inches

Conclusion:

Linear programming is a vital skill for anyone interested in Machine Learning/Data Science. In machine learning and deep learning, everything is about optimisation. At its very core you are trying to minimise some sort of loss function, and that itself is an optimisation problem. The main reason I suggest to learn linear programming is because it will help you understand the theoretical results of why convex/nonconvex optimisation methods word. ML algorithms are made up of convex or noncovex optimisation. The basic difference between the two categories is that in convex optimization there can be only one optimal solution, which is globally optimal or you might prove that there is no feasible solution to the problem, while in nonconvex optimization may have multiple locally optimal points and it can take a lot of time to identify whether the problem has no solution or if the solution is global. Hence, the efficiency in time of the convex optimization problem is much better. The main types of convex optimisation problems are logistic regression in a finite dimensional space or kernelised support vector machines in an infinite dimensional space. Nonconvex optimisation such as using deep neural networks is what most research is about nowadays. With this I hope you get a better understanding about Linear Programming and how to use it for other problems.