След края на Втората световна война Германия е разделена на четири окупационни зони: Американска, Британска, Френска и Съветска. Берлин също така е разделен на четири зони, но достъпът до него по суша и по вода минава изцяло през територията на Съветската окупационна зона.
На 24-ти юни 1948 Съветският Съюз блокира достъпа до Западен Берлин, с което започва една от най-значимите конфронтации по време на Студената война. Без достъп по суша или вода, западните съюзници започват да снабдяват града по въздух (Берлински въздушен мост).
За по-просто нека да приемем, че доставките до Берлин се извършват с два вида самолети: американски, които могат да поемат до 30 000 кубични фута товар (\approx 849 m^3), и британски самолети с капацитет до 20 000 кубични фута (\approx 566 m^3).
Колко британски и колко американски самолети да използва на ден въздушният мост, така че да достави до Берлин възможно най-голямо количество стоки?
\begin{align*} & x: \text{ брой американски самолети}\\ & y: \text{ брой британски самолети} \end{align*}
Общ товар, който x американски самолети и y британски самолети могат да доставят:
z(x, y) = 3 x + 2 y
\begin{align*} & x + y \leq 48 & \text{ (инфраструктура)} \\ & 16 x + 8 y \leq 512 & \text{ (персонал)} \\ & 9 x + 5 y \leq 300 & \text{ (бюджет)} \\ & x \geq 0 & \text{ (брой американски самолети)} \\ & y \geq 0 & \text{ (брой британски самолети)} \end{align*}
\begin{align*} & \max z(x, y) = 3 x + 2 y \\ & \text{при условията:} \\ & x + y \leq 48 \\ & 16 x + 8 y \leq 512 \\ & 9 x + 5 y \leq 300 \\ & x \geq 0 \\ & y \geq 0 \end{align*}
Преди да решим задачата, нека да я решим неформално в опростен вариант, като игнорираме второто и третото ограничение:
\begin{align*} \max 3 x + 2 y \\ x + y \leq 48 \\ x \geq 0 \\ y \geq 0 \end{align*}
Решение:
(x^* = 48, y^* = 0) \quad \text{и} \quad z^* = 3 \cdot 48 + 2 \cdot 0 = 144
За да решим задачата (и с трите ограничения) първо ще представим допустимото множество графично, като за целта ще начертаем правите към всяко от петте ограничения (включително ограниченията за неотрицателност). За да можем да начертаем правите са ни нужни по две точки от всяка права. Най-лесно можем да намерим пресечните им точки с двете оси (x и y).
За всички точки на оста y е изпълнено, че x = 0. Когато заместим с x = 0 в уравнението на всяка от правите ще намерим пресечните им точки с оста y. За всички точки от оста x важи, че y = 0, така че когато заместим с y = 0 в уравненията на правите ще получим пресечните им точки с оста x.
Права на първото ограничение (инфраструктура):
x + y = 48
За да намерим пресечните точки на правата (инфраструктура) с двете оси:
За да намерим пресечните точки на втората права (персонал) с двете оси:
16 x + 8 y = 512
По същия начин можем да намерим и две точки от правата на третото ограничение (бюджет): (0, 300 / 5) и (300 / 9, 0).
https://www.desmos.com/calculator/zvyvx248kr
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot([0, 48], [48, 0], label=r'$x + y = 48$ Инфрасруктура')
ax.plot([0, 32], [64, 0], label=r'$16 x + 8 y = 512$ Персонал')
ax.plot([0, 33.33], [60, 0], label=r'$9 x + 5 y = 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})")
plt.xlabel('x (американски самолети)')
plt.ylabel('y (британски самолети)')Text(0, 0.5, 'y (британски самолети)')
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^* = 15, y^* = 33). Това е възможно най-големият товар, който могат да пренесат самолетите при дадените ограничения. Този товар е равен на 30000 x^* + 20000 y^* = 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
x = m.addVar(vtype=GRB.INTEGER, lb = 0, name="американски самолети")
y = m.addVar(vtype=GRB.INTEGER, lb = 0, name="британски самолети")
# Целева функция
m.setObjective(3 * x + 2 * y, GRB.MAXIMIZE)
# Добавяне на ограниченията
m.addConstr(x + y <= 48, "Инфраструктура")
m.addConstr(16 * x + 8 * y <= 512, "Персонал")
m.addConstr(9 * x + 5 * y <= 300, "Бюджет")
# Неотрицателността на променливите е автоматично дефинирана в 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 |