Optimization Modelling in Python: Metaheuristics with constraints

Optimization modelling is one the most practical and widely used tools to find optimal or near-optimal solutions to complex decision-making problems.

In my previous post I gave example of very simple linear optimization problem with constraints, and provided exact solutions using several python libraries, namely PuLP, Pyomo, and SciPy.

In this article I will try to solve the same problem using new set of techniques called metaheuristics. In simple terms, metaheuristics is a generic, high-level procedure designed to find suboptimal solution, with little information or assumptions about given problem. Metaheuristics is gaining popularity in real world business problems and starting to be considered as real alternative to exact optimization methods. Here I will talk about two popular algorithms: Particle Swarm (PS) and Differential Evolution (DE).

Firstly, lets remind ourselves problem we are trying to solve. We have set of customers I = [1,2,3,4,5] and set of factories J = [1,2,3]. Each customer has some fixed product demand d_i and each factory has fixed production capacity M_j. We wont to distribute amount of goods sent from factories to customers in such a way that total transportation cost is minimal. Thus, this is integer minimization problem.

Problem has I binding constraints: sum of goods must be equal to customer demand. And J non-binding constraints: sum of goods sent must be no greater than factory production capacity.

Problem set up.

import time
import matplotlib.pyplot as plt
import numpy as np
d = {1:80 , 2:270 , 3:250 , 4:160 , 5:180}  # customer demand
M = {1:500 , 2:500 , 3:500}                 # factory capacity
I = [1,2,3,4,5]     # Customers
J = [1,2,3]         # Factories
cost = {(1,1):4,    (1,2):6,    (1,3):9,
(2,1):5,    (2,2):4,    (2,3):7,
(3,1):6,    (3,2):3,    (3,3):3,
(4,1):8,    (4,2):5,    (4,3):3,
(5,1):10,   (5,2):8,    (5,3):4
}              # transportation cost
# Transform cost dictionary into 2D array
cost2d = np.empty([len(I), len(J)])
for i in range(len(I)):
    for j in range(len(J)):
        cost2d[i,j] = cost[i+1,j+1]
# Variables bounds
n_vars = cost2d.size    # number of variables
bounds = 3*[(0,80), (0,270), (0,250), (0,160), (0,180)]

Exact solution

The exact solution of objective function is 3350 and nonzero variable values are as follows.

Exact solution using COIN-OR branch and cut solver

Converting constraints into penalties

As far as I know, being generic, metaheuristics does not support separate from objective function constraints definition. However, this is not a problem as we can convert constraints into penalties and include them into objective function.

For example, in given problem we have two sets of constraints: sum of goods must be equal to customer demand, and sum of goods must be no greater than factory capacity. Penalty in this case would be a unit of violation of constraint. For customer demand constraint, corresponding penalty could be written as max(0, abs(sum_j(goods_ij) - demand_i)**2. Because this constraint is binding, violation can be positive or negative, hence abs() in penalty. Raising penalty to power 2 is just to make it more severe.

For factory capacity constraint, corresponding penalty could be written as max(0, (sum_i(goods_ij) - capacity_j)**3. No abs() needed as constraint violation can be only positive in this case. Choice of power to raise penalty to is arbitrary. Power of 3 seemed to yield more consistently lower objective function values for this particular problem.

Basically, we are converting initial hard constraints into soft constraints, allowing for some of them being violated. Hence, sub-optimal solutions as a result.

Particle Swarm Optimization

In simple words PS is a population based method in which set of particles are exploring search space in stochastic fashion. Each particle movement has two components: local exploration - purely random local displacement, and social movement - displacement towards global best position of whole population. Thus, PS employs classical explore-exploit methodology.

The code snippet below is a vanilla PS algorithm from Dao, S.D. 2020.

if problem == 'minimize':
    initial_fitness = float("inf")
if problem == 'maximize':
    initial_fitness = -float("inf")
# Run Particle Swarm optimization
PSO(objective_function, bounds, n_particles, n_iterations)

Sub-optimal solution using Particle Swarm method

PS algorithm performs poorly for given problem. For a population of 200 particles with 500 iterations, for most of the runs value of objective function fluctuated around 5000. Just to remind you exact solution gives value 3350.

Even tuning explorative and social components of particle movement, or making them decay iteratively didn't improve result. PS method seems very prone to getting stuck in local minimas, from which particles cannot get out. Hence, quite large dispersion of objective function values.

In the PS snapshot above you can see penalty values for both constraints. Nonzero value for some penalties indicates violation of corresponding constraints. In this particular example, PS method allocated 260 and 154 units of goods to customers 2 and 4, while their demand was 270 and 160 correspondingly. Probably, not too bad.

Differential Evolution Optimization

Differential Evolution method is evolution based method that incorporates population based strategy. Similarly to PS, there are set of particles exploring search space stochastically. However, particle movement towards global best position is done differently. Instead of simply moving towards any local minima some ‘lucky’ particle may first fall into, DE employs particles mutation and crossover on each iteration. On each step, algorithm chooses one best particle, then creates new mutant particle from three randomly chosen ones, and finally mixes coordinates between local best particle and mutant. This so called gene/coordinates crossover is what prevents whole population from getting stuck in local minima.

bounds = 3*[(0,80), (0,270), (0,250), (0,160), (0,180)]
# Run Differential Evolution optimization
DE = list(de(objective_function, bounds, mut=0.8, crossp=0.7, popsize=200, its=500))
x = DE[-1][0]
ofunc_value = DE[-1][-1]
pen_cust_demand = \
    sum((max(0, abs(sum(x[idx: idx + len(J)]) - d[idx//len(J) +  1])))**2 for idx in range(0, cost2d.size, len(J)))
pen_fact_capacity = \
    sum((max(0, (sum(x[idx: idx + len(I)]) - M[idx//len(I) + 1])))**3 for idx in range(0, cost2d.size, len(I)))
print('RESULT:')
print('Objective function value:', ofunc_value)
print('Penalty customer demand:', pen_cust_demand)
print('Penalty factory capacity:', pen_fact_capacity)
print('Objective function value clean:', ofunc_value - pen_cust_demand - pen_fact_capacity)
EPS = 1.e-0
for idx,_ in enumerate(x):
    if x[idx] > EPS:
        print("sending quantity %10s from factory %3s to customer %3s" % (round(x[idx]), idx%len(J) + 1, idx//len(J) + 1))

Sub-optimal solution using Differential Evolution method

DE algorithm with 200 particles and 500 iterations has very good convergence and objective function value (3420 - 3450) is very close to exact solution (3350). Similarly to PS, some customers received little bit less goods than their demand.

Worth to note that both PS and DE may produce quite uneven variables values. I.e. sometimes some customers may be allocated 1 or 2 units of goods which is obviously not business meaningful.

In summary, metaheuristics is very promising set of optimization techniques that already can be used to real world problems. Similarly to machine learning models, these methods needs to be tuned to each particular problem. Given highly stochastic nature of these algorithms, each time they produce slightly different solutions.