За да разберем линейното оптимизиране, трябва да си припомним някои основни концепции от алгебрата и геометрията, които служат за построяването на нашите модели.
Основата на линейното оптимизиране са линейните зависимости. Най-често работим с уравнения на права в равнината (или хиперравнини в многомерното пространство).
Нека разгледаме пример за уравнение на права:
4x + 2y = -10
За да разберем геометричния смисъл на това уравнение, е удобно да го преобразуваме във вида y = mx + b, където m е наклонът (показва колко стръмна е правата и в каква посока върви. В случая, за всяка единица по оста X, стойността на y намалява с 2 единици. Отрицателният наклон означава, че правата "слиза" надолу.), а b е пресечната точка с оста Y (Това е точката, където правата пресича вертикалната ос Y. Когато x=0, y=-5).
Опростен вид: y = -2x - 5
Визуализация: Можете да експериментирате с графиките на функциите, използвайки Desmos Calculator. Опитайте да въведете горното уравнение, за да видите как промяната на коефициентите мести правата.
В задачите за оптимизация работим с ограничения, които дефинират лимитите на нашите ресурси (време, пари, суровини).
Ограниченията могат да бъдат:
Равенства: Изискват точно съвпадение. Изразяват строги технологични или логистични изисквания (напр. „Сумата от продуктите A и B трябва точно да запълни капацитета на контейнера”). Геометрично това са точките, лежащи точно върху правата (или равнината), което значително ограничава пространството от възможни решения в сравнение с неравенствата.
Неравенства: Дефинират диапазон (напр. "Разполагаме с 100 килограма от дадена суровина"). Геометрично те представляват полуравнини – всички точки от едната страна на правата.
Неотрицателност: В много практически задачи променливите описват количества ресурси, време, продукция или материали и поради физически или икономически смисъл не могат да приемат отрицателни стойности. Това обикновено ограничава решенията в първи квадрант, но е важно да се отбележи, че съществуват и неограничени променливи.
Допустимото множество е съвкупността от всички точки (решения), които удовлетворяват всички ограничения едновременно. В линейното оптимизиране това множество обикновено е многоъгълник (полигон). Ако допустимото множество е празно, задачата няма решение (невъзможно е да се удовлетворят всички условия). Оптималното решение при линейните програми винаги се намира във върхове на допустимото множество (може и да е повече от един връх).
Детерминираният модел приема, че всички параметри са известни със сигурност. Няма случайност в данните. За един и същ набор от входни данни, моделът винаги дава един и същ резултат (изход). Използват се, когато имаме фиксирано търсене, ясни разходи и точно определени ресурси. Обикновено са по-лесни за решаване изчислително.
Вероятностни / Стохастични модели включват елемент на несигурност. Поне един параметър се третира като случайна променлива и се описва чрез вероятностни разпределения (напр. нормално разпределение на търсенето). Изразяват се чрез очаквани стойности или вероятности, а не като едно твърдо число. Изискват по-сложен математически анализ и често симулации. Символичен пример: Марковските вериги моделират системи, които преминават от едно състояние в друго с определена вероятност - Визуализация на Марковски вериги
Задача 1: Програма за рационално хранене включва сено, смес и царевица. Тези храни са богати на белтъчини, калций и витамини. Количествата хранителни вещества съответно съдържащи се в 1 кг. от съответната храна са зададени в следната таблица:
| Белтъчини | Калций | Витамини | |
|---|---|---|---|
| Сено | 50 | 6 | 2 |
| Смес | 20 | 4 | 1 |
| Царевица | 180 | 3 | 1 |
Минимално необходимите количества белтъчини, калций и витамини (нормата на ден) са съответно 2000, 120 и 40. Съставете модел за определяне на оптималната комбинация от храни, при условие, че цените им са съответно 3, 2, 5.
\begin{align*} &\text{Неизвестни: } x_1, x_2, x_3 - \text{количества храна (в кг.)} \\ \\ &\text{Целева функция: } \\ &z(x) = 3x_1 + 2x_2 + 5x_3 \to \min \\ \\ &\text{Ограничения:} \\ &50x_1 + 20x_2 + 180x_3 \geq 2000 \\ &6x_1 + 4x_2 + 3x_3 \geq 120 \\ &2x_1 + x_2 + x_3 \geq 40 \\ &x_j \geq 0, \quad j = 1, 2, 3 \end{align*}
Задача 2: Предприятие произвежда 3 типа продукция – А, Б и В. За да произведе продукт А са необходими 2 минути машинно време, 30 минути труд и 7 €. допълнителни разходи. За продукт Б са необходими 2.5 минути машинно време, 40 минути труд и 10€. допълнителни разходи. За производството на продукт В са необходими 3 минути машинно време, 1 час труд и 13 €. допълнителни разходи. Продажната цена на продукт А е 12 €., на продукт Б – 14€., на продукт В – 20 €. Фирмата работи 5 дни в седмицата, с 2 смени по 7.5 часа дневно. Тя притежава 4 машини и за нея работят 50 човека на смяна. Седмичният бюджет за допълнителни разходи е 10 000 €. Фирмата цели максимизиране на приходите си от продажби. Съставете математическия модел.
\begin{align*} &\text{Целева функция:} \\ &z(x) = 12x_1 + 14x_2 + 20x_3 \to \max \\ \\ &\text{Ограничения за } X: \\ &2x_1 + 2.5x_2 + 3x_3 \leq 18\,000 \\ &30x_1 + 40x_2 + 60x_3 \leq 225\,000 \\ &7x_1 + 10x_2 + 13x_3 \leq 10\,000 \\ &x_j \geq 0, \quad j = 1, 2, 3 \end{align*}
Задача 3: Фирма - производител установява, че притежава излишък от производствени мощности. В резултат тя започва производство на 3 нови продукта. Разполагаемата мощност на машините, които се използват в производствения процес е 500, 300 и 150 часа съответно за машини А, Б, В. Броят часове, необходими за производството на всеки продукт е:
| Продукт 1 | Продукт 2 | Продукт 3 | |
|---|---|---|---|
| Машина А | 9 | 3 | 5 |
| Машина Б | 5 | 4 | 0 |
| Машина В | 3 | 0 | 2 |
Максималното количество продажби на продукт 3 не може да надвиши 20 единици. Цените на продуктите са съответно 30, 12, 15.
\begin{align*} &\text{Целева функция:} \\ &z(x) = 30x_1 + 12x_2 + 15x_3 \to \max \\ \\ &\text{Ограничения за } X: \\ &9x_1 + 3x_2 + 5x_3 \leq 500 \\ &5x_1 + 4x_2 \leq 300 \\ &3x_1 + 2x_3 \leq 150 \\ &x_3 \leq 20 \\ &x_j \geq 0, \quad j = 1, 2, 3 \end{align*}
След края на Втората световна война Германия е разделена на четири окупационни зони: Американска, Британска, Френска и Съветска. Берлин също така е разделен на четири зони, но достъпът до него по суша и по вода минава изцяло през територията на Съветската окупационна зона.
На 24-ти юни 1948 Съветският Съюз блокира достъпа до Западен Берлин, с което започва една от най-значимите конфронтации по време на Студената война. Без достъп по суша или вода, западните съюзници започват да снабдяват града по въздух (Берлински въздушен мост).
За по-просто нека да приемем, че доставките до Берлин се извършват с два вида самолети: американски, които могат да поемат до 30 000 кубични фута товар (\approx 849 m^3), и британски самолети с капацитет до 20 000 кубични фута (\approx 566 m^3).
Поради ограничения в инфраструктурата на ден могат да летят най-много 48 самолета (независимо от кой вид).
За всеки полет на американски самолет има нужда от 16 души персонал, двойно повече от броя нужен за британските самолети. Общо на разположение са 512 души на ден.
Разходите за гориво и поддръжка на самолетите възлизат на 9000 долара за полет на американски самолет и на 5000 долара за полет на британски самолет. Поради бюджетни ограничения общите разходи не могат да надхвърлят 300 000 долара.
Колко британски и колко американски самолети да използва на ден въздушният мост, така че да достави до Берлин възможно най-голямо количество стоки?
\begin{align*} & x_A: \text{ брой американски самолети}\\ & x_B: \text{ брой британски самолети} \end{align*}
Общ товар, който x_A американски самолети и x_B британски самолети могат да доставят:
z(x_A, x_B) = 3 x_A + 2 x_B
\begin{align*} & x_A + x_B \leq 48 & \text{ (инфраструктура)} \\ & 16 x_A + 8 x_B \leq 512 & \text{ (персонал)} \\ & 9 x_A + 5 x_B \leq 300 & \text{ (бюджет)} \\ & x_A \geq 0 & \text{ (брой американски самолети)} \\ & x_B \geq 0 & \text{ (брой британски самолети)} \end{align*}
\begin{align*} & \max z(x_A, x_B) = 3 x_A + 2 x_B \\ & \text{при условията:} \\ & x_A + x_B \leq 48 \\ & 16 x_A + 8 x_B \leq 512 \\ & 9 x_A + 5 x_B \leq 300 \\ & x_A \geq 0 \\ & x_B \geq 0 \end{align*}
За да решим задачата (и с трите ограничения) първо ще представим допустимото множество графично, като за целта ще начертаем правите към всяко от петте ограничения (включително ограниченията за неотрицателност). За да можем да начертаем правите са ни нужни по две точки от всяка права. Най-лесно можем да намерим пресечните им точки с двете оси (x_A и x_B).
За всички точки на оста x_B е изпълнено, че x_A = 0. Когато заместим с x_A = 0 в уравнението на всяка от правите ще намерим пресечните им точки с оста x_B. За всички точки от оста x_A важи, че x_B = 0, така че когато заместим с x_B = 0 в уравненията на правите ще получим пресечните им точки с оста x_A.
Права на първото ограничение (инфраструктура):
x_A + x_B = 48
За да намерим пресечните точки на правата (инфраструктура) с двете оси:
За да намерим пресечните точки на втората права (персонал) с двете оси:
16 x_A + 8 x_B = 512
По същия начин можем да намерим и две точки от правата на третото ограничение (бюджет): (0, 300 / 5) и (300 / 9, 0).
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot([0, 48], [48, 0], label=r'$x_A + x_B = 48$ Инфрасруктура')
ax.plot([0, 32], [64, 0], label=r'$16 x_A + 8 x_B = 512$ Персонал')
ax.plot([0, 33.33], [60, 0], label=r'$9 x_A + 5 x_B = 300$ Бюджет')
ax.legend(loc=0)
import polytope as pc
import numpy as np
A = np.array([[1, 1], [16, 8], [9, 5], [-1, 0], [0, -1]])
b = np.array([48, 512, 300, 0, 0])
P = pc.Polytope(A, b)
P_extr = pc.extreme(P)
P_extr`polytope` failed to import `cvxopt.glpk`.
will use `scipy.optimize.linprog`array([[32., -0.],
[20., 24.],
[15., 33.],
[-0., 48.],
[-0., -0.]])Допустимото множество се състои от всички точки в полигона (0, 0), (32, 0), (20, 24), (15, 33), (0, 48). Изчислението на пресечните точки (20, 24) и (15, 33).
from matplotlib import pyplot as plt
plt.plot(P_extr[:, 0], P_extr[:, 1], '-o')
for i, p in enumerate(P_extr):
plt.text(p[0], p[1], f"({p[0]:0.1f}, {p[1]:0.1f})")
def z(xA, xB):
return 3 * xA + 2 * xB
# (0, 0)
print("z(0, 0) = ", z(0, 0))
# (32, 0)
print("z(32, 0) = ", z(32, 0))
# (20, 24)
print("z(20, 24) = ", z(20, 24))
# (15, 33)
print("z(15, 33) = ", z(15, 33))
# (0, 48)
print("z(0, 48) = ", z(0, 48))z(0, 0) = 0
z(32, 0) = 96
z(20, 24) = 108
z(15, 33) = 111
z(0, 48) = 96Оптималната комбинация от американски и британски самолети е (x^*_A = 15, x^*_B = 33). Това е възможно най-големият товар, който могат да пренесат самолетите при дадените ограничения. Този товар е равен на $30000 x^*_{A} + 20000 x^*_{B} = 30000 \cdot 15 + 20000 \cdot 33 = 1 110 000 кубични фута, използвайки 15 американски и 33 британски самолета.
В момента решихме задачата, като изчислихме върховете на допустимото множество и пресметнахме целевата функция във всеки от тях. Решението на максимизационната задача беше върхът с най-висока стойност на целевата функция.
За съжаление този подход е приложим само за много малки задачи и няма практическа стойност. Причината за това е, че броят на върховете на допустимото множество нараства много бързо с увеличаване на броя на ограниченията и променливите. Горна граница за броя на върховете на допустимото множество е дадена от биномния коефициент:
\binom{n}{k} = \frac{n!}{k!(n-k)!}
където n е броят на променливите, а k е броят на ограниченията. Можем да пресметнем горната граница за броя на върховете на допустимото за различни комбинации на n и k:
from scipy.special import comb
print(f"n = 20, k = 15, Максимален брой върхове = {comb(20, 15):,.0f}")
print(f"n = 40, k = 25, Максимален брой върхове = {comb(40, 25):,.0f}")
print(f"n = 50, k = 25, Максимален брой върхове = {comb(50, 25):,.0f}")n = 20, k = 15, Максимален брой върхове = 15,504
n = 40, k = 25, Максимален брой върхове = 40,225,345,056
n = 50, k = 25, Максимален брой върхове = 126,410,606,437,752Изчисляването на целевата функция за всички върхове е непостижимо дори за модерни компютри. Ако приемем, че компютърът може да обработи 1 милиард върха на секунда, за 80 променливи и 40 ограничения ще му трябва много време:
comb_n = comb(80, 45)
years = comb_n / (1e9 * 60 * 60 * 24 * 365)
print(f"{years:,.0f} години")1,836,017 годиниimport gurobipy as gp
from gurobipy import GRB
import pandas as pd
# Съставяне на нов модел
m = gp.Model("Berlin Airlift")
m.Params.LogToConsole = 0
xA = m.addVar(vtype=GRB.INTEGER, lb = 0, name="американски самолети")
xB = m.addVar(vtype=GRB.INTEGER, lb = 0, name="британски самолети")
# Целева функция
m.setObjective(3 * xA + 2 * xB, GRB.MAXIMIZE)
# Добавяне на ограниченията
m.addConstr(xA + xB <= 48, "Infrastructure")
m.addConstr(16 * xA + 8 * xB <= 512, "Staff")
m.addConstr(9 * xA + 5 * xB <= 300, "Budget")
# Неотрицателността на променливите е автоматично дефинирана в m.addVar
# Решаване на модела
m.optimize()
print("Максимален обем: ", 1e4 * m.objVal, "кубични фута")
print("Използвайки комбинация от:")
# Отпечатване на резултата в таблица
df = pd.DataFrame(columns=["Променлива", "Стойност"])
for v in m.getVars():
df.loc[len(df)] = [v.varName, v.x]
dfRestricted license - for non-production use only - expires 2027-11-29
Set parameter LogToConsole to value 0
Максимален обем: 1110000.0 кубични фута
Използвайки комбинация от:| Променлива | Стойност | |
|---|---|---|
| 0 | американски самолети | 15.0 |
| 1 | британски самолети | 33.0 |