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 massdefobjective(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 runsol[:,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)**2returnobj

## 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.