Optimization deals with selecting the best option among a number of possible choices that are feasible or don't violate constraints. Python can be used to optimize parameters in a model to best fit data, increase profitability of a potential engineering design, or meet some other type of objective that can be described mathematically with variables and equations.

Mathematical optimization problems may include equality constraints (e.g. =), inequality constraints (e.g. <, <=, >, >=), objective functions, algebraic equations, differential equations, continuous variables, discrete or integer variables, etc. One example of an optimization problem from a benchmark test set is the Hock Schittkowski problem #71.

$$\min x_1 x_4 \left(x_1 + x_2 + x_3\right) + x_3$$ $$\mathrm{s.t.} \quad x_1 x_2 x_3 x_4 \ge 25$$ $$x_1^2 + x_2^2 + x_3^2 + x_4^2 = 40$$ $$1\le x_1, x_2, x_3, x_4 \le 5$$ $$x_0 = (1,5,5,1)$$

This problem has a nonlinear objective that the optimizer attempts to minimize. The variable values at the optimal solution are subject to (s.t.) both equality (=40) and inequality (>25) constraints. The product of the four variables must be greater than 25 while the sum of squares of the variables must also equal 40. In addition, all variables must be between 1 and 5 and the initial guess is x_{1} = 1, x_{2} = 5, x_{3} = 5, and x_{4} = 1.

import numpy as np

from scipy.optimize import minimize

def objective(x):

return x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]

def constraint1(x):

return x[0]*x[1]*x[2]*x[3]-25.0

def constraint2(x):

sum_eq = 40.0

for i in range(4):

sum_eq = sum_eq - x[i]**2

return sum_eq

# initial guesses

n = 4

x0 = np.zeros(n)

x0[0] = 1.0

x0[1] = 5.0

x0[2] = 5.0

x0[3] = 1.0

# show initial objective

print('Initial Objective: ' + str(objective(x0)))

# optimize

b = (1.0,5.0)

bnds = (b, b, b, b)

con1 = {'type': 'ineq', 'fun': constraint1}

con2 = {'type': 'eq', 'fun': constraint2}

cons = ([con1,con2])

solution = minimize(objective,x0,method='SLSQP',\

bounds=bnds,constraints=cons)

x = solution.x

# show final objective

print('Final Objective: ' + str(objective(x)))

# print solution

print('Solution')

print('x1 = ' + str(x[0]))

print('x2 = ' + str(x[1]))

print('x3 = ' + str(x[2]))

print('x4 = ' + str(x[3]))

from scipy.optimize import minimize

def objective(x):

return x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]

def constraint1(x):

return x[0]*x[1]*x[2]*x[3]-25.0

def constraint2(x):

sum_eq = 40.0

for i in range(4):

sum_eq = sum_eq - x[i]**2

return sum_eq

# initial guesses

n = 4

x0 = np.zeros(n)

x0[0] = 1.0

x0[1] = 5.0

x0[2] = 5.0

x0[3] = 1.0

# show initial objective

print('Initial Objective: ' + str(objective(x0)))

# optimize

b = (1.0,5.0)

bnds = (b, b, b, b)

con1 = {'type': 'ineq', 'fun': constraint1}

con2 = {'type': 'eq', 'fun': constraint2}

cons = ([con1,con2])

solution = minimize(objective,x0,method='SLSQP',\

bounds=bnds,constraints=cons)

x = solution.x

# show final objective

print('Final Objective: ' + str(objective(x)))

# print solution

print('Solution')

print('x1 = ' + str(x[0]))

print('x2 = ' + str(x[1]))

print('x3 = ' + str(x[2]))

print('x4 = ' + str(x[3]))

from gekko import GEKKO

import numpy as np

#Initialize Model

m = GEKKO()

#help(m)

#define parameter

eq = m.Param(value=40)

#initialize variables

x1,x2,x3,x4 = [m.Var() for i in range(4)]

#initial values

x1.value = 1

x2.value = 5

x3.value = 5

x4.value = 1

# lower bounds

x1.lower = 1

x2.lower = 1

x3.lower = 1

x4.lower = 1

# upper bounds

x1.upper = 5

x2.upper = 5

x3.upper = 5

x4.upper = 5

#Equations

m.Equation(x1*x2*x3*x4>=25)

m.Equation(x1**2+x2**2+x3**2+x4**2==eq)

#Objective

m.Obj(x1*x4*(x1+x2+x3)+x3)

#Set global options

m.options.IMODE = 3 #steady state optimization

#Solve simulation

m.solve()

#Results

print('')

print('Results')

print('x1: ' + str(x1.value))

print('x2: ' + str(x2.value))

print('x3: ' + str(x3.value))

print('x4: ' + str(x4.value))

import numpy as np

#Initialize Model

m = GEKKO()

#help(m)

#define parameter

eq = m.Param(value=40)

#initialize variables

x1,x2,x3,x4 = [m.Var() for i in range(4)]

#initial values

x1.value = 1

x2.value = 5

x3.value = 5

x4.value = 1

# lower bounds

x1.lower = 1

x2.lower = 1

x3.lower = 1

x4.lower = 1

# upper bounds

x1.upper = 5

x2.upper = 5

x3.upper = 5

x4.upper = 5

#Equations

m.Equation(x1*x2*x3*x4>=25)

m.Equation(x1**2+x2**2+x3**2+x4**2==eq)

#Objective

m.Obj(x1*x4*(x1+x2+x3)+x3)

#Set global options

m.options.IMODE = 3 #steady state optimization

#Solve simulation

m.solve()

#Results

print('')

print('Results')

print('x1: ' + str(x1.value))

print('x2: ' + str(x2.value))

print('x3: ' + str(x3.value))

print('x4: ' + str(x4.value))

This tutorial can also be completed with nonlinear programming optimizers that are available with the Excel Solver and MATLAB Optimization Toolbox. Click on the appropriate link for additional information and source code.