Решение задач линейного программирования с использованием Python

Зачем решать экстремальные задачи

На практике очень часто возникают задачи, для решения которых используются методы оптимизации. В обычной жизни при множественном выборе, например, подарков к новому годы мы интуитивно решаем задачу минимальных затрат при заданном качестве покупок.

К сожалению, не всегда можно положиться на интуицию. Допустим Вы сотрудник коммерческой фирмы и отвечаете за рекламу. Затраты на рекламу в месяц не должны превышать 10 000 денежных единиц (д.е). Минута радиорекламы стоит 5 д.е., а телерекламы 90 д.е. Фирма намерена использовать радиорекламу в два раза чаще чем телерекламу. Практика показывает, что 1 минута телерекламы обеспечивает объём продаж в 30 раз больший чем 1 минута радиорекламы.

Перед Вами стоит задача определить такое распределение средств между двумя упомянутыми видами рекламы при котором объём продаж фирмы будет максимальным. Вы сначала выберите переменные, а именно месячный объём в минутах на телерекламу — x1, а на радиорекламу --x2. Теперь не трудно составить следующую систему:

30x1+x2 –увеличение продаж от рекламы;
90x1+5x2 <=10 000 – ограничение средств;
x2=3x1 – соотношение времён радио и теле рекламы.

Теперь на время забудем о рекламе и постараемся обобщить изложенное. Таких задач, как приваленная, может быть много, но они имеют следующие общие признаки. Обязательным есть наличие линейно зависящей от переменных функции цели, в нашем случае это 30x1+x2, которая при найденных значениях входящих переменных должна иметь единственное максимальное значение. При этом условие не отрицательности входящих переменных выполняется автоматически. Далее следуют опять-таки линейные равенства и неравенства в количестве, зависящем от наличия условий. Вот мы и сформулировали одну группу задач линейного программирования.

Другую большую группу задач линейного программирования, рассмотрим на примере так называемой транспортной задачи. Допустим Вы сотрудник коммерческой фирмы, которая оказывает транспортные услуги. Есть поставщики товара со складами в разных трёх городах, причём объёмы однородной продукции на этих складах соответственно равны a1, a2, a3. Есть и потребители в других трёх городах которым нужно привести товар от поставщиков в объёмах b1, b2, b3 соответственно. Известны также стоимости доставки с1÷с9 товаров от поставщиков к потребителям, согласно таблице.

Если обозначить через x1…xn количество перевозимого груза, тогда функцией цели будет общая стоимость перевозки:

F(x)=c1*x1+c2*x2+c3*x3+c4*x4+c5*x5+c6*x6+c7*x7+c8*x8+c9*x9.

Условия, которые записываться. в виде неравенств:

x1+x2+x3<=20 – больше чем есть у поставщика не возьмёшь
x4+x5+x6<=45
x7+x8+x9<=30

Условия, которые записываться. в виде равенств:

x1+x4+x7=b1– сколько надо столько и привезём
x2+x5+x8=b2
x3+x6+x9=b3

Тут дополнительно нужны условия не отрицательности переменных x поскольку они по смыслу не отрицательны и ищется минимум F(x). Эти неравенства не приводим.

Теперь Вы знаете как строить функции цели и условия для основных задач линейного программирования. Но когда Вы прочитаете в специальной литературе про геометрический, симплекс, искусственного базиса методы решения указанных задач Вы бросили и рекламу и логистику. Но ведь можно найти простое и понятное решение на Python.

Выбор библиотек Python для решения типовых задач линейного программирования

Для линейного программирования в Python мне известны три специализированные библиотеки. Рассмотрим решение обеих приведенных задач на каждой из библиотек. Кроме интерфейса и результатов оптимизации будем оценивать и быстродействие. Поскольку нам нужно только качественное отличие в быстродействии воспользуемся для этого самым простым листингом усредняя результаты каждого запуска программ.

Оптимизация с библиотекой pulp [1].

from pulp import *
import time
start = time.time()
x1 = pulp.LpVariable("x1", lowBound=0)
x2 = pulp.LpVariable("x2", lowBound=0)
problem = pulp.LpProblem('0',pulp.LpMaximize)
problem += 30*x1 +x2, "Функция цели"
problem += 90*x1+ 5*x2 <= 10000, "1"
problem +=x2 ==3*x1, "2"
problem.solve()
print ("Результат:")
for variable in problem.variables():
    print (variable.name, "=", variable.varValue)
print ("Прибыль:")
print (value(problem.objective))
stop = time.time()
print ("Время :")
print(stop - start)

В лис тенге программы уже знакомые нам соотношения для максимальной прибыли от рекламы 30*x1+x2, условия ограничения затрат, помеченные для сравнения «1». Мы не забыли и об отношении времён использования радио и теле рекламы, помеченные в лис тенге как «2». Назначение других операторов очевидны, Подробности можно прочесть в [1].

Результаты решения задачи оптимизации с использованием pulp.

Результат:
x1 = 95.238095
x2 = 285.71429
Прибыль:
3142.85714
Время:
0.10001182556152344

В итоге мы получили значения времён рекламы, при которых ожидаемая прибыль от её использования будет максимальна.

Оптимизация с библиотекой cvxopt [2].

from cvxopt.modeling import variable, op
import time
start = time.time()
x = variable(2, 'x')
z=-(30*x[0] +1*x[1])#Функция цели
mass1 = (90*x[0] + 5*x[1]  <= 10000) #"1"
mass2 = (3*x[0] -x[1] == 0) # "2"
x_non_negative = (x >= 0) #"3"    
problem =op(z,[mass1,mass2,x_non_negative])
problem.solve(solver='glpk')  
problem.status
print ("Прибыль:")
print(abs(problem.objective.value()[0]))
print ("Результат:")
print(x.value)
stop = time.time()
print ("Время :")
print(stop - start)

По структуре программа аналогична предыдущей, но имеются два существенных отличия. Во-первых, библиотека cvxopt настроена на поиск минимума функции цели, а не на максимум. Поэтому целевая функция взята с отрицательным знаком минус -(30*x[0] +1*x[1]). Полученное вследствие этого отрицательное её значение выведено по абсолютной величине. Во-вторых, введено ограничение на не отрицательность переменных- non_negative. Повлияло ли это на результат мы сейчас у видим.

Результаты решения задачи оптимизации с использованием cvxopt.

Прибыль:
3142.857142857143
Результат:
[ 9.52e+01]
[ 2.86e+02]
Время:
0.041656494140625

Никаких существенных изменений в сравнении с библиотекой pulp не произошло за исключением время работы программы.

Оптимизация с библиотекой scipy. optimize [3].

from scipy.optimize import linprog
import time
start = time.time()
c = [-30,-1] #Функция цели
A_ub = [[90,5]]  #'1'   
b_ub = [10000]#'1'   
A_eq = [[3,-1]] #'2'   
b_eq = [0] #'2'   
print (linprog(c, A_ub, b_ub, A_eq, b_eq))
stop = time.time()
print ("Время :")
print(stop - start)

Достаточно беглого взгляда на листинг, чтобы понять, что мы имеем дело с принципиально иным подходом к вводу данных. Хотя приведенные в листингах цифры помогают прояснить принцип организации данных, путём сравнения, всё же приведу пояснения. Список c = [-30,-1] содержит коэффициенты функции цели с обратным знаком, поскольку linprog () ищет минимум. Матрица A_ub содержит коэффициенты при переменных для условий в виде неравенств. Для нашей задачи это 90x1+5x2 <=10000. Значения в правой части неравнства-1000, помещается в список b_ub. Матрица A_eq содержит коэффициенты при переменных для условий в виде равенств. Для нашей задачи 3x1-x2=0, причём ноль в правой части, помещается в список b_eq.

Результаты решения задачи оптимизации с использованием scipy. optimize.

Fun: -3142.8571428571431
message: 'Optimization terminated successfully.'
nit: 2
slack: array([ 0.])
status: 0
success: True
x: array ([ 95.23809524, 285.71428571])

Время:
0.03020191192626953

Здесь более подробно выведены результаты расчётов, но сами результаты те же что и в предыдущих библиотеках.

По результатам решения задачи «О рекламе» можно сделать промежуточный вывод, о том что использование библиотеки scipy. optimize обеспечивает большее быстродействие и рациональную форму исходных данных. Однако без результатов решения транспортной задачи окончательный вывод делать рано.

Привожу решение транспортной задачи, но уже без подробных пояснений, поскольку основные этапы решения уже подробно описаны.

Оптимизация с библиотекой pulp.

from pulp import *
import time
start = time.time()
x1 = pulp.LpVariable("x1", lowBound=0)
x2 = pulp.LpVariable("x2", lowBound=0)
x3 = pulp.LpVariable("x3", lowBound=0)
x4 = pulp.LpVariable("x4", lowBound=0)
x5 = pulp.LpVariable("x5", lowBound=0)
x6 = pulp.LpVariable("x6", lowBound=0)
x7 = pulp.LpVariable("x7", lowBound=0)
x8 = pulp.LpVariable("x8", lowBound=0)
x9 = pulp.LpVariable("x9", lowBound=0)
problem = pulp.LpProblem('0',pulp.LpMaximize)
problem += -7*x1 - 3*x2 - 6* x3 - 4*x4 - 8*x5 -2* x6-1*x7- 5*x8-9* x9, "Функция цели"
problem +=x1 + x2 +x3<= 74,"1" 
problem +=x4 + x5 +x6 <= 40, "2"
problem +=x7 + x8+ x9 <= 36, "3"
problem +=x1+ x4+ x7 == 20, "4"
problem +=x2+x5+ x8 == 45, "5"
problem +=x3 + x6+x9 == 30, "6"                     
problem.solve()
print ("Результат:")
for variable in problem.variables():
    print (variable.name, "=", variable.varValue)
print ("Стоимость доставки:")
print (abs(value(problem.objective)))
stop = time.time()
print ("Время :")
print(stop - start)

Решение транспортной задачи требует минимизации затрат по доставке, поэтому функция цели вводиться со знаком минус, а выводиться по абсолютной величине.

Результаты решения транспортной задачи с использованием pulp.

Результат:
x1 = 0.0
x2 = 45.0
x3 = 0.0
x4 = 0.0
x5 = 0.0
x6 = 30.0
x7 = 20.0
x8 = 0.0
x9 = 0.0
Стоимость доставки:
215.0
Время:
0.19006609916687012

Оптимизация с библиотекой cvxopt.

from cvxopt.modeling import variable, op
import time
start = time.time()
x = variable(9, 'x')
z=(7*x[0] + 3*x[1] +6* x[2] +4*x[3] + 8*x[4] +2* x[5]+x[6] + 5*x[7] +9* x[8])
mass1 = (x[0] + x[1] +x[2] <= 74)
mass2 = (x[3] + x[4] +x[5] <= 40)
mass3 = (x[6] + x[7] + x[8] <= 36)
mass4 = (x[0] + x[3] + x[6] == 20)
mass5 = (x[1] +x[4] + x[7] == 45)
mass6 = (x[2] + x[5] + x[8] == 30)
x_non_negative = (x >= 0)    
problem =op(z,[mass1,mass2,mass3,mass4 ,mass5,mass6, x_non_negative])
problem.solve(solver='glpk')  
problem.status
print("Результат:")
print(x.value)
print("Стоимость доставки:")
print(problem.objective.value()[0])
stop = time.time()
print ("Время :")
print(stop - start)

Результаты решения транспортной задачи с использованием cvxopt.

Результат:
[ 0.00e+00]
[ 4.50e+01]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 3.00e+01]
[ 2.00e+01]
[ 0.00e+00]
[ 0.00e+00]
Стоимость доставки:
215.0
Время :
0.03001546859741211

Оптимизация с библиотекой scipy. optimize.

from scipy.optimize import linprog	
import time
start = time.time()
c = [7, 3, 6,4,8,2,1,5,9]
A_ub = [[1,1,1,0,0,0,0,0,0],
               [0,0,0,1,1,1,0,0,0],
               [0,0,0,0,0,0,1,1,1]] 
b_ub = [74,40,36] 
A_eq = [[1,0,0,1,0,0,1,0,0],
               [0,1,0,0,1,0,0,1,0],
               [0,0,1,0,0,1,0,0,1]] 
b_eq = [20,45,30] 
print(linprog(c, A_ub, b_ub, A_eq, b_eq))
stop = time.time()
print ("Время :")
print(stop - start)

Результаты решения транспортной задачи с использованием scipy optimize.

fun: 215.0
message: 'Optimization terminated successfully.'
nit: 9
slack: array([ 29., 10., 16.])
status: 0
success: True
x: array([ 0., 45., 0., 0., 0., 30., 20., 0., 0.])
Время:
0.009982585906982422

Анализ решения двух типовых задач линейного программирования с помощью трёх библиотек аналогичного назначения не вызывает сомнения в выборе библиотеки scipy. optimize, как лидера по компактности ввода данных и быстродействию.

Что нового для использования библиотеки scipy. optimize при решении задач линейного программирования

Получение из исходной матрицы, списка для функции цели, а также заполнение матриц неравенств A_ub и равенств A_eq программно, позволит упростить работу по вводу данных и увеличить размерность исходной матрицы. Это позволить решать реальные задачи оптимизации. Рассмотрим, как это можно сделать на простом демонстрационном примере, не претендующем на идеальность кода.

import numpy as np
from scipy.optimize import linprog
b_ub = [74,40,36] 
b_eq = [20,45,30] 
A=np.array([[7, 3,6],[4,8,2],[1,5,9]])
m, n = A.shape
c=list(np.reshape(A,n*m))# Преобразование матрицы A в список c.
A_ub= np.zeros([m,m*n])
for i in np.arange(0,m,1):# Заполнение матрицы условий –неравенств.
         for j in np.arange(0,n*m,1):
                  if i*n<=j<=n+i*n-1:
                        A_ub  [i,j]=1
A_eq= np.zeros([m,m*n])
for i in np.arange(0,m,1):# Заполнение матрицы условий –равенств.
         k=0
         for j in np.arange(0,n*m,1):
                  if j==k*n+i:
                           A_eq [i,j]=1
                           k=k+1
print(linprog(c, A_ub, b_ub, A_eq, b_eq))

Теперь вводиться только сама матрица A и списки правых частей b_ub неравенств и b_ub – равенств.

Результат рaботы программы предсказуем.
fun: 215.0
message: 'Optimization terminated successfully.'
nit: 9
slack: array([ 29., 10., 16.])
status: 0
success: True
x: array([ 0., 45., 0., 0., 0., 30., 20., 0., 0.])

Вывод частный

При решении задач линейного программирования целесообразно использовать библиотеку scipy.optimize по методике, приведенной в статье, а матрицы условий заполнять програмно. При этом Вам не понадобятся специальные знания о методах решения оптимизационных задач.

Вывод общий

В последнее время появились разные библиотеки Python решающие аналогичные задачи. Решение о применении той или иной библиотеки часто носит субъективный характер. Поэтому целесообразно проводить их сравнительный анализ для области решаемых Вами задач.

Ссылки

1. pythonhosted.org/PuLP
2. cvxopt.org/userguide/modeling.html
3. docs.scipy.org/doc/scipy/reference/tutorial/optimize.html

Автор: Scorobey

Источник