Зачем решать экстремальные задачи
На практике очень часто возникают задачи, для решения которых используются методы оптимизации. В обычной жизни при множественном выборе, например, подарков к новому году мы интуитивно решаем задачу минимальных затрат при заданном качестве покупок.
К сожалению, не всегда можно положиться на интуицию. Допустим Вы сотрудник коммерческой фирмы и отвечаете за рекламу. Затраты на рекламу в месяц не должны превышать 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 решающие аналогичные задачи. Решение о применении той или иной библиотеки часто носит субъективный характер. Поэтому целесообразно проводить их сравнительный анализ для области решаемых Вами задач.
Ссылки
Поделиться с друзьями
Комментарии (8)
Dark_Daiver
11.06.2017 10:01+2Статью пробежал по диагонали, но очень бросилось в глаза то, что тестовая задача слишком маленькая для тестирования производительности библиотек. Тест который выполняется сотые доли секунды это совсем не показательно для интерпретируемого языка.
Я совершенно не специалист в оптимизациях, но все же я не уверен что корректно сравнивать scipy linprog и cvxopt.
Первый — просто вызов симплекс солвера по готовым матрицам.
https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.optimize.linprog.html
cvxopt же — целый DSL (т.е. вам не нужно в ручную подготавливать данные для солвера) позволяющий решать намного больший класс задач.
Так же, в примере с cvxopt использовался солвер glpk (https://www.gnu.org/software/glpk/), который вроде как предназначен для работы с задачами большой размерности («The GLPK (GNU Linear Programming Kit) package is intended for solving large-scale linear programming (LP)»). Т.е. в статье сравнивается теплое с мягким.
Scorobey
12.06.2017 03:42Цель статьи не тестирование библиотек, а сравнение интерфейсных возможностей решателей для двух конкретных групп задач линейного программирования.Различные формы представления исходных данных с использованием матрицы для формирования функции цели и без влияет на аппаратное время выполнения программы.При усреднении времен такое сравнение корректно и оправдано.В статье не рассматривались возможности библиотек в целом, поэтому утверждение о том что cvxopt позволяет решать намного больший класс задач отношение к статье не имеет. Кроме того само утверждения голословно — какой класс, каких задач.
Спасибо за комментарий.Dark_Daiver
12.06.2017 08:27Окей, давайте по порядку.
>Цель статьи не тестирование библиотек, а сравнение интерфейсных возможностей решателей для двух конкретных групп задач линейного программирования
В таком случае, на это стоило делать акцент.
К тому же, мне кажется что если мы сравниваем библиотеки по способу задания задачи, cvxopt/pulp явно предпочтительней. Почему? Потому что они позволяют писать значительно более читаемый и поддерживаемый код.
Ограничения добавляются/удаляются в cvxopt одной строчкой. Функционал так же очень легко модифицируется. Фактически вы просто пишете формулы описывающие ограничения и функционал.
Если же матрицы сформированы вручную, то это разумеется более эффективно, но читать и работать с таким кодом довольно неприятно. Как минимум вам придется писать комментарии
Собственно код из вашей статьи.
Это cvxopt. Если нормально назвать переменные (x1 — moneyForTVAd x2 — moneyForRadioAd) и ф-ии (z — negativeSalesIncr и т.д.), то понять какую именно задачу вы решаете, сможет даже человек не знакомый с предметной областью
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"
Это та же самая задача, но решаемая при помощи linprog (если я не ошибаюсь, конечно).
Работать можно и с этим, но согласитесь, это значительно менее приятно чем работать с кодом выше.
c = [-30,-1] #Функция цели A_ub = [[90,5]] #'1' b_ub = [10000]#'1' A_eq = [[3,-1]] #'2' b_eq = [0] #'2'
Идем дальше
>Различные формы представления исходных данных с использованием матрицы для формирования функции цели и без влияет на аппаратное время выполнения программы.
Да, это разумеется так.
>При усреднении времен такое сравнение корректно и оправдано.
И это тоже верно, вот только я не увидел чтобы вы решали запускали решение LP задач тысячу раз, чтобы говорить об усреднении. Но это не самое страшное. Проблема в том, что на задачах очень малой размерности на время выполнения начинает сильно влиять различный .Overhead. Overhead от вызова самих питоновских функций (если анализ задачи оптимизации в cvxopt выполняется на питоне, то он может занимать больше времени чем сама оптимизация), overhead от вызова C/C++ кода, и т.д.
Более того, как я и писал выше, linprog и cvxopt (в вашем случае) используют разные алгоритмы оптимизации заточенные под разные случаи. Симплекс-метод хорошо подходит для небольших задач (как ваши), GLPK же (если я правильно понял) использует итеративный алгоритм и заточен под большие задачи (т.е. сотни тысяч переменных/ограничений). Подозреваю, что у cvxopt есть и другие солверы, однако у вас вызывается именно этот. Собственно на этом можно остановиться, сравнение производительности в такой ситуации смысла не имеет.
>Кроме того само утверждения голословно — какой класс, каких задач.
Выпуклые оптимизации же =) CVXOPT — ConVeX OPTimization.
Merlen_Gross
Python 2.* и вырвиглазное форматирование — классика.
Andy_U
Python 3.x — print не оператор, а функция. Ну, пробел перед скобкой… Строчки «from __future__ import print_function» нет.