Распределение ресурсов в условиях неопределенности

Управление доходами (англ. Revenue management, сокращённо RM) звучит, как что-то очень скучное. Максимизация прибыли, усиление конкурентоспособности, эффективное планирование и бюджетирование, улучшение принятия решений, устойчивое развитие. Разве не скука? Также всё это управление доходами может показаться циничным, ведь в таких сферах, как медицина и образование, это зачастую становится причиной несправедливых решений.

Однако! Благодаря RM компании развиваются. Развитие компаний — это развитие всего рынка. Развитие рынка — это рост экономики. Рост экономики — это увеличение: налоговых поступлений, количества рабочих мест, качества жизни и благополучия общества.

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

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

RM как дисциплина достаточно универсальна. Абстрактно RM можно представить как механизм, который перерабатывает данные в оптимальные решения. Часто RM воспринимается как машинное обучение (ML) в чистом виде. Некоторые компании даже отправляют своих сотрудников на курсы по ML, желая на самом деле обучить их RM, и это, несомненно, приводит к казусам, например, ошибочное восприятие предиктивных моделей как конечного результата. Машинное обучение и датасаенс (Data science, наука о данных, сокращённо DS) в целом очень важны, но являются лишь первым главным этапом решения задач RM. Представим, что данные — это нефть. Их нужно переработать в керосин, на котором будет работать двигатель математического программирования и исcледования операций. Чем лучше качество топлива и тюнинг двигателя — тем лучше решения, получаемые на выходе.

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

Задача оптимизации распределения ресурсов в условиях неопределенности состоит в том, что мы оптимизируем математическое ожидание возможных значений целевой функции, а не какое-то ее конкретное значение. Концептуально это является единственным отличием от детерминированных задач. Тем не менее, это отличие влечет за собой много сложностей и неочевидных нюансов. Рассмотрим базовый принцип оптимизации на нескольких простых примерах.

Оптимальный объем одного типа ресурса

Предположим, что некоторая авиакомпания решила расширить свою маршрутную сеть по некоторому направлению с помощью кодшеринга и уже определилась с оперирующим перевозчиком. Осталось только понять — на каких условиях заключать соглашение и сколько мест нужно выкупить у оперирующего перевозчика.

В первую очередь авиакомпании необходимо выяснить параметры спроса по новому направлению. Байесовский вывод является, пожалуй, самым разумным способом сделать это. Даже если нет никаких данных, то благодаря байесовской статистике можно начинать с относительно разумных предположений о параметрах спроса, а далее уточнять их по мере появления данных о перевозках. Пока предположим, что авиакомпания решила задачу вероятностного программирования (не путать со стохастическим программированием) и в ее распоряжении есть PMF (Probability Mass Function) распределения спроса.

Распределение спроса может аппроксимироваться какими-то известными законами распределения или быть эмпирическим, имея сложную форму. Главное — это наличие PMF. Предположим, авиакомпания выяснила, что распределение спроса наилучшим образом аппроксимируется бета-биномиальным распределением $\varphi(x) = f(x \; | \; n, a, b, \mathrm{loc})$ , в котором , , и $\mathrm{loc} = 10$ :

Python

# Выполним необходимые импорты

import numpy as np
import pandas as pd
import networkx as nx
np.random.seed(42)
from scipy import stats
from scipy.stats import uniform, betabinom, randint, bernoulli, binom, poisson
import matplotlib
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 6, 3.5
rcParams['figure.dpi'] = 140
%config InlineBackend.figure_format = 'png'
import seaborn as sns
#sns.set()
#from matplotlib import cm
#sns.set_style("whitegrid")

n, a, b, loc = 50, 4, 3, 10

x = np.arange(10, 60)
y = betabinom.pmf(x, n, a, b, loc)

plt.plot(x, y, 'bo', ms=5)
plt.vlines(x, 0, y, colors='b', lw=1, alpha=0.5)
plt.axvline(np.sum(x*y), color='r', label = 'mean')
plt.legend()
plt.xlabel('Quantity demanded')
plt.ylabel('pmf')
plt.title('Distribution of demand')
plt.ylim(0, 0.041);

На графике видно, что распределение является несимметричным.

Возможные убытки, скорее всего, также будут несимметричными. Например:

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

Обозначим убытки в первом случае $k_{1}$ , а во втором случае $k_{2}$ . Очевидно, что авиакомпания покупает места у оперирующего перевозчика дешевле, чем продает, следовательно $k_{1} < k_{2}$ . Пусть данные виды убытков будут иметь следующие значения:

$k_{1}$ = 1000 у.е.
$k_{2}$ = 1150 у.е.

Python

# k_1 - убытки от одного пустого кресла:
k_1 = 1000

# k_2 - убытки от одного нехватившего кресла:
k_2 = 1150

Теперь авиакомпании необходимо определиться с количеством мест , которое ей необходимо выкупить у оперирующего перевозчика. Имея информацию о распределении спроса и возможных убытках, крайне трудно принять верное решение, поэтому действовать приходится "на глаз". Предположим, что в авиакомпании решили выбирать между математическим ожиданием, медианой и модой. Допустим, выбор пал на значение = 40, близкое к матожиданию:

Python

# c - выделенное количество мест:
c = 40

Учитывая возможные значения спроса, необходимо определить возможные значения количества пустых и недостающих мест $X_{u}$ и $X_{m}$ соответственно:

Python

# количества незанятых мест:
X_u = c - x[x <= c]
# количества недостающих мест:
X_m = x[x >= c] - c

f, ax1 = plt.subplots()

# Распределение количества свободных мест:
ax1.bar(x[x <= c], X_u, color='g', alpha=0.5, label='Unoccupied')
# Распределение количества недостающих мест:
ax1.bar(x[x >= c], X_m, color='r', alpha=0.5, label='Missing')
ax1.legend()
ax1.set_xlabel('Demand')
ax1.set_ylabel('Number of seats') 
ax1.set_title('Possible numbers of missing and unoccupied seats')
ax1.grid(False)

ax2 = ax1.twinx()

ax2.set_ylim(0, max(y) + 0.005)
ax2.plot(x, y, 'bo', ms=5)
ax2.vlines(x, 0, y, colors='b', lw=1, alpha=0.5)
ax2.set_ylabel('pdf')
ax2.grid(False);

Поскольку спрос имеет вероятностную природу, то убытки также будут случайными. Хорошее представление о распределении убытков можно получить, если взглянуть на взвешенные значения $X_{u}$ и $X_{m}$ :

Python

f, ax1 = plt.subplots()

# Распределение количества незанятых мест:
U_dist = X_u * y[x <= c]
ax1.bar(x[x <= c], U_dist, color='g', alpha=0.5, label='Unoccupied')
U_pdf = U_dist / np.sum(U_dist)
E_u = np.sum(x[x <= c] * U_pdf)
ax1.axvline(E_u, color='green')
# Распределение количества необходимых мест:
M_dist = X_m * y[x >= c]
ax1.bar(x[x >= c], M_dist, color='r', alpha=0.5, label='Missing')
ax1.legend()
M_pdf = M_dist / np.sum(M_dist)
E_m = np.sum(x[x >= c] * M_pdf)
ax1.axvline(E_m, color='red')
ax1.set_xlabel('Demand')
ax1.set_ylabel('Number of seats') 
ax1.set_title(r'Weighted values of $X_{u}$ and $X_{m}$')
ax1.grid(False)

ax2 = ax1.twinx()



ax2.set_ylim(0, max(y) + 0.005)
ax2.plot(x, y, 'bo', ms=5)
ax2.vlines(x, 0, y, colors='b', lw=1, alpha=0.5)
ax2.set_ylabel('pdf')
ax2.grid(False);

Сразу заметно, что распределение значений несимметрично. Может показаться, что в постановке задачи есть некоторое заблуждение: ведь непроданные кресла являются реальными убытками, а потенциальные клиенты, которым не хватило места — это потенциальная недополученная прибыль. Однако и в первом, и во втором случае речь идет о деньгах, которые теряет авиакомпания. И чем меньше эти потери, тем лучше. Таким образом, убытки от пустых и недостающих кресел $L_{u}$ и $L_{m}$ можно рассматривать как некий возможный убыток, обусловленный двумя причинами и имеющий собственное распределение:

Python

# убыток от незанятых мест:
L_u = k_1 * (c - x[x <= c])
# убыток от недостающих мест:
L_m = k_2 * (x[x >= c] - c)

L = np.hstack((L_u, L_m[1:]))

plt.bar(x, L*y, color='tab:blue')
plt.xlabel('Demand')
plt.ylabel('Loss share')
plt.title('Distribution of shares from the average loss');

Средний общий убыток при таком количестве купленных мест составит:

$\mathbb {E}(L) = \sum^{c}_{i=\mathrm{loc}}x^{(u)}_{i}\varphi(x^{(u)}_{i}) + \sum^{n}_{i=c+1}x^{(m)}_{i}\varphi(x^{(m)}_{i}),$

где $x^{(u)}_{i}$ и $x^{(m)}_{i}$ — конкретные значения количества пустых и недостающих кресел соответственно. Для выбранного авиакомпанией значения = 40 значение $\mathbb {E}(L)$ составит 8167 у.е.

Python

np.sum(L*y)

Результат

8167.534771308343

Вычислим все возможные значения средних общих убытков для разных значений :

Python

Loss = []
for c in x:
    L_u = k_1 * (c - x[x <= c])
    L_m = k_2 * (x[x >= c] - c)
    L = np.hstack((L_u, L_m[1:]))
    Loss.append(np.sum(L * y))

plt.plot(x, Loss)
plt.plot(x, Loss, 'bo', ms=5)
plt.vlines(x, min(Loss) - 1000, Loss, colors='b', lw=1, alpha=0.5)
plt.plot(x[np.argmin(Loss)], min(Loss), 'o', ms=5, color='red')
plt.vlines(x[np.argmin(Loss)], min(Loss) - 1000, min(Loss), colors='red', lw=2, alpha=0.5)

plt.xlabel(r'$c$')
plt.ylabel('loss (c.u.)')
plt.ylim(min(Loss) - 1000, 34000)
plt.title('Dependence of average total losses on c');

Взглянув на результаты на графике, можно подумать, что авиакомпания угадала с выбранным значением = 40, так как $\mathrm{argmin}(L(c))$ также равен 40. Однако минимизация средних общих убытков вовсе не означает максимизацию средней прибыли. Прибыль в данном случае равна разности среднего дохода и среднего общего убытка.

$\mathbb {E} \, ( \, \mathrm {Profit} \, ) = \mathbb {E} \, ( \, \mathrm {Income} \, ) - \mathbb {E} \, ( \, \mathrm {Loss} \, )$

Средний доход $\mathbb {E} (I)$ ограничивается либо величиной спроса, либо количеством мест, которые авиакомпания купит у оперирующего перевозчика. Чтобы найти среднее количество занаятых кресел, необходимо принять во внимание, что места могут быть использованы либо полностью, либо частично:

$\mathbb {E} (X) = \sum^{c}_{i=\mathrm{loc}}x_{i}\varphi(x_{i}) + \sum^{n}_{i=c+1}c\varphi(x_{i}).$

Например, для = 40 среднее количество занятых мест $\mathbb {E} (X)$ составит около 35.5 штук.

Python

# Теоретический расчет E(X):
p_up = 40*(1 - betabinom.cdf(39, n, a, b, loc))
p_up

p_down = sum(y[x < 40] * x[x < 40])
p_down

p_up + p_down

Результат

35.42930863685582

Python

# Экспериментальное значение E(X):
D = betabinom.rvs(n, a, b, loc, size=10000)
D[D > 40] = 40
D.mean()

Результат

35.4861

С учетом, что прибыль от проданного места равняется $k_{2}$ , мы теперь можем вычислить средний доход $\mathbb {E} (I)$ при различных значениях :

$\mathbb {E} (I) = k_{2} \left ( \sum^{c}_{i=\mathrm{loc}}x_{i}\varphi(x_{i}) + \sum^{n}_{i=c+1}c\varphi(x_{i}) \right ).$

Python

e_profit = []

for c in x:
    p_up = c*(1 - betabinom.cdf(c - 1, n, a, b, loc))
    p_down = sum(y[x < c] * x[x < c])
    p = p_up + p_down
    e_profit.append(p * k_2)
    
plt.plot(x, e_profit, 'bo-')
plt.vlines(x, min(e_profit), e_profit, colors='b', lw=1, alpha=0.5)
plt.ylim(min(e_profit)-500, max(e_profit)+3000)
plt.xlabel(r'$c$')
plt.ylabel('income (c.u.)')
plt.title(r'Dependence of average income on $c$');

Построим зависимость средней прибыли $\mathbb {E} \, ( \, \mathrm {Profit} \, )$ от :

Python

Loss = np.array(Loss)
e_profit = np.array(e_profit)
Z = e_profit - Loss

plt.plot(x, Z)
plt.plot(x, Z, 'bo', ms=5)
plt.vlines(x, min(Z) - 1000, Z, colors='b', lw=1, alpha=0.5)
plt.plot(x[np.argmax(Z)], max(Z), 'o', ms=5, color='red')
plt.vlines(x[np.argmax(Z)], min(Z) - 1000, max(Z),
           label=r'$c$=44', colors='red', lw=2, alpha=0.5)
plt.ylim(min(Z) - 500, max(Z) + 3000)
plt.legend()
plt.xlabel(r'$c$')
plt.ylabel('profit (c.u.)')
plt.title(r'Dependence of average profit on $c$');

Максимальное значение средней прибыли наблюдается при = 44 и составляет 33622 у.е.

Python

np.max(Z)

Результат

33622.00462011786

Значение средней общей прибыли при изначальном значении = 40 составит 32576 у.е.

Python

Z[40 - 10]

Результат

32576.17016107585

В процентном выражении разность составляет 3.2%:

Python

100 * (Z[44 - 10] - Z[40 - 10]) / Z[40 - 10]

Результат

3.210427910557898

Получение 3% дополнительной прибыли — это весьма отлично. Будем честны: дополнительная прибыль — это ведь всегда отлично, даже если она составляет всего 0.5%. Однако интересно, что если бы авиакомпания не использовала оптимизацию вовсе, то ее решение можно было бы интерпретировать как случайное угадывание. Парадокс заключается в том, что принятое авиакомпанией решение может быть даже как-то обосновано и восприниматься абсолютно рациональным. Ведь неудачу всегда можно списать на случайную природу обстоятельств, а успех присвоить себе.

Если предположить, что значения выбираются случайным образом из интервала (который может показаться наиболее разумным), то процент потерянной средней прибыли будет выглядеть следующим образом:

Python

lp = [100 * (Z[44 - 10] - Z[i - 10]) / Z[i - 10] for i in range(35, 51)]

plt.plot(np.r_[35:51], lp, 'bo-', zorder=10)
plt.vlines(np.r_[35:51], 0, lp, colors='b', lw=1, alpha=0.5)
plt.xlabel(r'$c$')
plt.ylabel('%', rotation = 0, labelpad=10)
plt.title('Loss of average total profit')
plt.axhline(0, color='0.5');

Если не разбираться в сути проблемы, то потеря прибыли может принимать большие значения. Самое печальное, что об этом очень долго можно даже не подозревать.

Важность моделирования

Источником всех сложностей управления ресурсами является асимметрия: типов убытков, распределения спроса, средних доходов и средних общих убытков. Любая асимметрия — это источник нелинейности, из-за которой аналитический вывод общего вида целевой функции и ее ограничений может быть как затруднен, так ошибочен. Человеческий фактор в данном случае можно исключить только посредством моделирования.

Много ошибок было допущено уже в самом начале решения даже этой простой задачи. Дело вовсе не в уровне образования или навыков решающего задачу, а в том, что есть множество нюансов, которые трудно заметить без эксперимента. Ошибки являются неотъемлемой частью исследовательского процесса. Казалось бы, это прописная истина, но только задумайтесь, к чему может привести сумма человеческого фактора и токсичного отношения к ошибкам. Как часто такие факторы суммируются на практике? А теперь представьте, что на практике чаще всего встречается что-то вроде:

"Заткнись и считай так, как положено!"
"Этим должно заниматься целое НИИ, а не какой-то один умник!"
"Марьиванна лучше знает, ведь она все 45 лет так делает, и ничего!"

Смех смехом, но это абсолютно реальные жизненные ситуация. И это я еще даже не коснулся темы оборудования! Кто знает, что такое вероятностное программирование, тот знает, сколько надо "жечь железа". Добавитть сюда стохастическое программирование и умножить необходимые вычислительные мощности еще на два.

Важность моделирования может показаться сомнительной, потому что пример действительно очень простой. Давайте задумаемся: у оперирующего перевозчика может быть две кабины с разными классами обслуживания, и авиакомпания могла бы захотеть купить места в каждой из них. Также авиакомпания могла бы создать несколько классов бронирования с группами тарифов в каждом из них, а значит убытки также пришлось бы считывать более сложным образом. По каждому классу бронирования пришлось бы выводить отдельную PMF распределения спроса на билеты данного подкласса. Количество объектов в задаче и количество связей между ними значительно увеличится, а с учетом того, что все эти связи, как правило, нелинейны, то и ошибку допустить крайне легко.

Моделирование — это не просто хорошая практика, порой это единственный способ добиться оптимизации.

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

Оптимальные объемы нескольких типов ресурсов

Теперь рассмотрим, как определять оптимальный объем для нескольких типов ресурсов.

Предположим, что некоторой авиакомпании необходимо разделить места в самолете, вмещающего 100 кресел, на два класса обслуживания: бизнес и эконом. Авиакомпания уже обладает необходимыми данными и получила на их основе соответствующие PMF распределения спроса, которые выглядят следующим образом:

Python

# Параметры распределения спроса
# на места в кабине бизнес класса:
n_biz, a_biz, b_biz, loc_biz = 35, 5, 6, 0

# Параметры распределения спроса
# на места в кабине эконом класса:
n_eco, a_eco, b_eco, loc_eco = 120, 6, 5, 0

# pmf бизнес класса:
x_biz = np.arange(0, 35)
pmf_biz = betabinom.pmf(x_biz, n_biz, a_biz, b_biz, loc_biz)

# pmf эконом класса:
x_eco = np.arange(0, 120)
pmf_eco = betabinom.pmf(x_eco, n_eco, a_eco, b_eco, loc_eco)

f, ax = plt.subplots(1, 2, figsize=(12, 4))

ax[0].plot(x_biz, pmf_biz, 'bo', ms=3)
ax[0].vlines(x_biz, 0, pmf_biz, colors='b', lw=1, alpha=0.5)
ax[0].set_title('pmf of business class demand')

ax[1].plot(x_eco, pmf_eco, 'bo', ms=3)
ax[1].vlines(x_eco, 0, pmf_eco, colors='b', lw=1, alpha=0.5)
ax[1].set_title('pmf of demand for economy class');

В качестве примера предположим, что авиакомпания уже определилась со значениями количества кресел в каждом классе (или уже использовала их ранее):

$c^{(biz)}$ = 16;
$c^{(eco)}$ = 84.

Python

# N - общая вместимость самолета:
N = 100

# c_biz - выделенное количество мест для бизнесс класса:
c_biz = 16

# c_eco - выделенное количество мест для эконом класса:
c_eco = 84

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

Python

# количества незанятых мест в бизнесс классе:
X_biz_u = c_biz - x_biz[x_biz <= c_biz]
# количества необходимых мест в бизнесс классе:
X_biz_m = x_biz[x_biz >= c_biz] - c_biz

# количества незанятых мест в эконом классе:
X_eco_u = c_eco - x_eco[x_eco <= c_eco]
# количества необходимых мест в эконом классе:
X_eco_m = x_eco[x_eco >= c_eco] - c_eco


f, (ax1, ax3)= plt.subplots(2, 1, figsize=(11, 8))

# Распределение количества незанятых мест в бизнесс классе:
U_biz_dist = X_biz_u * pmf_biz[x_biz <= c_biz]
ax1.bar(x_biz[x_biz <= c_biz], U_biz_dist, color='g', alpha=0.5, label='Unoccupied')
U_biz_pdf = U_biz_dist / np.sum(U_biz_dist)
E_biz_u = np.sum(x_biz[x_biz <= c_biz] * U_biz_pdf)
ax1.axvline(E_biz_u, color='green')

# Распределение количества необходимых мест в бизнесс классе:
M_biz_dist = X_biz_m * pmf_biz[x_biz >= c_biz]
ax1.bar(x_biz[x_biz >= c_biz], M_biz_dist, color='r', alpha=0.5, label='Missing')
M_biz_pdf = M_biz_dist / np.sum(M_biz_dist)
E_biz_m = np.sum(x_biz[x_biz >= c_biz] * M_biz_pdf)
ax1.axvline(E_biz_m, color='red')

ax1.set_ylabel('Number of seats') 
ax1.grid(False)
ax1.set_title('Business Class')

ax2 = ax1.twinx()

ax2.set_ylim(0, max(pmf_biz) + 0.005)
ax2.plot(x_biz, pmf_biz, 'bo', ms=5)
ax2.vlines(x_biz, 0, pmf_biz, colors='b', lw=1, alpha=0.5)
ax2.set_ylabel('pdf')
ax2.grid(False)


# Распределение количества незанятых мест в эконом классе:
U_eco_dist = X_eco_u * pmf_eco[x_eco <= c_eco]
ax3.bar(x_eco[x_eco <= c_eco], U_eco_dist, color='g', alpha=0.5, label='Unoccupied')
U_eco_pdf = U_eco_dist / np.sum(U_eco_dist)
E_eco_u = np.sum(x_eco[x_eco <= c_eco] * U_eco_pdf)
ax3.axvline(E_eco_u, color='green')

# Распределение количества необходимых мест в эконом классе:
M_eco_dist = X_eco_m * pmf_eco[x_eco >= c_eco]
ax3.bar(x_eco[x_eco >= c_eco], M_eco_dist, color='r', alpha=0.5, label='Missing')
M_eco_pdf = M_eco_dist / np.sum(M_eco_dist)
E_eco_m = np.sum(x_eco[x_eco >= c_eco] * M_eco_pdf)
ax3.axvline(E_eco_m, color='red')

ax3.set_xlabel('Demand')
ax3.set_ylabel('Number of seats') 
ax3.grid(False)
ax3.set_title('Economy class')

ax4 = ax3.twinx()

ax4.set_ylim(0, max(pmf_eco) + 0.005)
ax4.plot(x_eco, pmf_eco, 'bo', ms=5)
ax4.vlines(x_eco, 0, pmf_eco, colors='b', lw=1, alpha=0.5)
ax4.set_ylabel('pdf')
ax4.grid(False);

Однако теперь емкость самолета является фиксированной, и повлиять на убытки от недостающих кресел никак не удастся. Иными словами, теперь задача максимизации прибыли будет совпадать с задачей минимизации убытков от недостающих кресел.

Пусть $k^{\mathrm {(class)}}_{u}$ и $k^{\mathrm {(class)}}_{m}$ — это убытки от нераспроданных и недостающих кресел соответственно. Пусть они имеют следующие значения:

$k^{\mathrm {(biz)}}_{u}$ = 2000 у.е.;
$k^{\mathrm {(biz)}}_{m}$ = 2300 у.е.;
$k^{\mathrm {(eco)}}_{u}$ = 1000 у.е.;
$k^{\mathrm {(eco)}}_{m}$ = 1100 у.е..

Python

# k_biz_u - убытки от одного пустого кресла в бизнесс классе:
k_biz_u = 2000    # 2300

# k_biz_m - убытки от одного нехватившего кресла в бизнесс классе:
k_biz_m = 2300    # 2000

# k_eco_u - убытки от одного пустого кресла в эконом классе:
k_eco_u = 1000    # 1100

# k_eco_m - убытки от одного нехватившего кресла в эконом классе:
k_eco_m = 1100    # 1000

Для вычисления средней общей прибыли понадобятся следующие величины:

$\mathbb {E} (I)$ — средние общие доходы;
$\mathbb {E}(L_{u})$ — средние общие убытки от нераспроданных кресел;
$\mathbb {E}(L_{m})$ — средние общие убытки от недостающих кресел.

От распределения мест по классам эти величины будут зависеть следующим образом:

Python

def sales_param(c_biz, c_eco, len_d):
    d_biz = betabinom.rvs(n_biz, a_biz, b_biz, loc_biz, size = len_d)
    diff_biz = d_biz - c_biz
    
    if len(diff_biz[diff_biz < 0]) == 0:
        eu_biz = 0
    else:
        eu_biz = np.mean(np.abs(diff_biz[diff_biz < 0]) * k_biz_u)
        
    if len(diff_biz[diff_biz > 0]) == 0:
        em_biz = 0
    else:
        em_biz = np.mean(diff_biz[diff_biz > 0] * k_biz_m)
        
    d_biz[d_biz > c_biz] = c_biz
    p_biz = np.mean(d_biz * k_biz_m)

    
    d_eco = betabinom.rvs(n_eco, a_eco, b_eco, loc_eco, size = len_d)
    diff_eco = d_eco - c_eco
    
    if len(diff_eco[diff_eco < 0]) == 0:
        eu_eco = 0
    else:
        eu_eco = np.mean(np.abs(diff_eco[diff_eco < 0]) * k_eco_u)
        
    if len(diff_eco[diff_eco > 0]) == 0:
        em_eco = 0
    else:
        em_eco = np.mean(diff_eco[diff_eco > 0] * k_eco_m)
    d_eco[d_eco > c_eco] = c_eco
    p_eco = np.mean(d_eco * k_eco_m)

    eu = eu_biz + eu_eco
    em = em_biz + em_eco
    p = p_biz + p_eco

    return p, eu, em

params = np.array([sales_param(i, 100 - i, 500000) for i in range(0, 35)])
P = params[:, 0]
Eu = params[:, 1]
Em = params[:, 2]

plt.plot(P, marker='o', label=r'$\mathbb {E} (I)$')
plt.plot(Eu, marker='o', label=r'$\mathbb {E}(L_{u})$')
plt.plot(Em, marker='o', label=r'$\mathbb {E}(L_{m})$')
plt.xlabel(r'$c^{biz}$ ($c^{eco} = N - c^{biz}$)')
plt.ylabel('c.u.')
plt.legend()
plt.title('Dependence of the average total income and losses\n\
on the distribution of seats');

При моделировании можно заметить, что решения для максимизации прибыли, которая рассчитывается как $\mathbb {E} (I) - \mathbb {E}(L_{u})$ , и решения для максимизации прибыли, которая рассчитывается как $\mathbb {E} (I) - \mathbb {E}(L_{u}) - \mathbb {E}(L_{m})$ , могут отличаться:

Python

plt.plot(P-Eu, label=r'$\mathbb {E} (I) - \mathbb {E}(L_{u})$')
plt.plot(P-Eu-Em, label=r'$\mathbb {E} (I) - \mathbb {E}(L_{u}) - \mathbb {E}(L_{m})$')
plt.plot(np.argmax(P-Eu), max(P-Eu), 'ro')
plt.axvline(np.argmax(P-Eu), color='0.7', ls='--')
plt.plot(np.argmax(P-Eu-Em), max(P-Eu-Em), 'ro')
plt.axvline(np.argmax(P-Eu-Em), color='0.7', ls='--')
plt.legend()
plt.xlabel(r'$c^{biz}$ ($c^{eco} = N - c^{biz}$)')
plt.ylabel('profit (c.u.)')
plt.title('Optimal solutions for different types of profit');

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

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

Обозначим возможные значения спроса на места в бизнес и эконом классах как $d^{\mathrm {biz}}_{i}$ и $d^{\mathrm {eco}}_{j}$ , тогда возможные общие убытки от пустых кресел будут равны:

$L_{ij} = k^{\mathrm {biz}}_{u}(с^{\mathrm {biz}} - d^{\mathrm {biz}}_{i}) + k^{\mathrm {eco}}_{u}(с^{\mathrm {eco}} - d^{\mathrm {eco}}_{j}).$

А их распределение будет выглядеть следующим образом:

Python

# Распределение убытков в бизнес классе:
L_biz_dist = k_biz_u * (c_biz - x_biz[x_biz <= c_biz])

# Распределение убытков в эконом классе:
L_eco_dist = k_eco_u * (c_eco - x_eco[x_eco <= c_eco])

L_biz_dist = L_biz_dist.reshape(L_biz_dist.size, -1)
L_eco_dist = L_eco_dist.reshape(1, L_eco_dist.size)
L_joint_dist = L_biz_dist + L_eco_dist

im = plt.imshow(L_joint_dist)
plt.xlabel(r'$d^{\mathrm {eco}}_{j}$')
plt.ylabel(r'$d^{\mathrm {biz}}_{i}$', rotation=0)
plt.title('Distribution of possible total losses')

plt.colorbar(im, orientation='horizontal');

Совместная PMF спроса может быть вычислена как:

$\varphi(d^{\mathrm {(biz)}}_{i}, d^{\mathrm {(eco)}}_{j}) = \frac{\varphi^{\mathrm {(biz)}}(d^{\mathrm {(biz)}}_{i}) \cdot \varphi^{\mathrm {(eco)}}(d^{\mathrm {(eco)}}_{j})}{\sum^{c^{\mathrm {(biz)}}}_{i=loc^{\mathrm {(biz)}}}\sum^{c^{\mathrm {(eco)}}}_{i=loc^{\mathrm {(eco)}}}\varphi^{\mathrm {(biz)}}(d^{\mathrm {(biz)}}_{i}) \cdot \varphi^{\mathrm {(eco)}}(d^{\mathrm {(eco)}}_{j})}$

Python

pmf_biz_2d = pmf_biz[x_biz <= c_biz].reshape(pmf_biz[x_biz <= c_biz].size, -1)
pmf_eco_2d = pmf_eco[x_eco <= c_eco].reshape(1, pmf_eco[x_eco <= c_eco].size)
pmf_joint_dist = pmf_biz_2d * pmf_eco_2d
pmf_joint_dist = pmf_joint_dist/ np.sum(pmf_joint_dist)

im = plt.imshow(pmf_joint_dist)
plt.xlabel(r'$d^{\mathrm {eco}}_{j}$')
plt.ylabel(r'$d^{\mathrm {biz}}_{i}$', rotation=0)
plt.title('Joint distribution of demand')

plt.colorbar(im, orientation='horizontal');

Взвешенные общие убытки примут следующий вид:

Python

plt.imshow(L_joint_dist * pmf_joint_dist);

А само значение среднего общего убытка при выбранной комбинации мест в каждом классе будет следующим:

$\mathbb {E} (L(c^{\mathrm {(biz)}},c^{\mathrm {(eco)}})) = \sum^{c^{\mathrm {(biz)}}}_{i=loc^{\mathrm {(biz)}}}\sum^{c^{\mathrm {(eco)}}}_{i=loc^{\mathrm {(eco)}}}L_{ij} \cdot \varphi(d^{\mathrm {(biz)}}_{i}, d^{\mathrm {(eco)}}_{j}).$

Например, для выбранных авиакомпанией значений $c^{(biz)}$ = 16 и $c^{(eco)}$ = 84 средний общий убыток будет равен 32274 у.е.

Python

E_loss = np.sum(L_joint_dist * pmf_joint_dist)
E_loss

Результат

32273.514081521025

Создадим функцию для вычисления средних общих убытков:

Python

def loss(c_biz, c_eco, k_biz_u, k_biz_m, k_eco_u, k_eco_m):
    L_biz_dist = k_biz_u * (c_biz - x_biz[x_biz <= c_biz])

    # Распределение убытков в эконом классе:
    L_eco_dist = k_eco_u * (c_eco - x_eco[x_eco <= c_eco])

    L_biz_dist = L_biz_dist.reshape(L_biz_dist.size, -1)
    L_eco_dist = L_eco_dist.reshape(1, L_eco_dist.size)
    L_joint_dist = L_biz_dist + L_eco_dist
    
    # Совместная плотность вероятности:
    pmf_biz_2d = pmf_biz[x_biz <= c_biz].reshape(pmf_biz[x_biz <= c_biz].size, -1)
    pmf_eco_2d = pmf_eco[x_eco <= c_eco].reshape(1, pmf_eco[x_eco <= c_eco].size)
    pmf_joint_dist = pmf_biz_2d * pmf_eco_2d
    pmf_joint_dist = pmf_joint_dist/ np.sum(pmf_joint_dist)
    
    E_loss = np.sum(L_joint_dist * pmf_joint_dist)
    
    return E_loss

Средний общий доход вычисляется по следующей формуле:

$\mathbb {E} (I(c^{\mathrm {(biz)}},c^{\mathrm {(eco)}})) = k^{\mathrm {(biz)}} \left ( \sum^{c^{\mathrm {(biz)}}}_{i=loc^{\mathrm {(biz)}}} d^{\mathrm {(biz)}} \varphi^{\mathrm {(biz)}}(d^{\mathrm {(biz)}}_{i}) + \sum^{n^{\mathrm {(biz)}}}_{i=c^{\mathrm {(biz)}} + 1} c^{\mathrm {(biz)}} \cdot \varphi^{\mathrm {(biz)}}(d^{\mathrm {(biz)}}_{i}) \right ) + \\ + \; k^{\mathrm {(eco)}} \left ( \sum^{c^{\mathrm {(eco)}}}_{i=loc^{\mathrm {(eco)}}} d^{\mathrm {(eco)}} \varphi^{\mathrm {(eco)}}(d^{\mathrm {(eco)}}_{i}) + \sum^{n^{\mathrm {(eco)}}}_{i=c^{\mathrm {(eco)}} + 1} c^{\mathrm {(eco)}} \cdot \varphi^{\mathrm {(eco)}}(d^{\mathrm {(eco)}}_{i}) \right ).$

Для выбранного авиакомпанией распределения количества кресел по классам он будет равен 101846 у.е.

Python

p_up = c_biz * (1 - betabinom.cdf(c_biz - 1, n_biz, a_biz, b_biz, loc_biz))
p_down = sum(pmf_biz[x_biz < c_biz] * x_biz[x_biz < c_biz])
p_biz = k_biz_m*(p_up + p_down)

p_up = c_eco * (1 - betabinom.cdf(c_eco - 1, n_eco, a_eco, b_eco, loc_eco))
p_down = sum(pmf_eco[x_eco < c_eco] * x_eco[x_eco < c_eco])
p_eco = k_eco_m*(p_up + p_down)
p_biz + p_eco

Результат

101845.81112969195

Значение средней прибыли будет равно:

$\mathbb {E} (\mathrm {Profit}(c^{\mathrm {(biz)}},c^{\mathrm {(eco)}})) = \mathbb {E} (I(c^{\mathrm {(biz)}},c^{\mathrm {(eco)}})) - \mathbb {E} (L(c^{\mathrm {(biz)}},c^{\mathrm {(eco)}})).$

Необходимо найти значения $c^{(biz)}$ и $c^{(eco)}$ , при которых $\mathbb {E} (\mathrm {Profit})$ достигнет максимума при условии, что $c^{(biz)} + c^{(eco)} = 100$ .

Python

NP = np.full((101, 101), -np.inf)

for c_biz in x_biz:
    c_eco = 100 - c_biz
    p_up = c_biz * (1 - betabinom.cdf(c_biz - 1, n_biz, a_biz, b_biz, loc_biz))
    p_down = sum(pmf_biz[x_biz < c_biz] * x_biz[x_biz < c_biz])
    p_biz = k_biz_m*(p_up + p_down)
    p_up = c_eco * (1 - betabinom.cdf(c_eco - 1, n_eco, a_eco, b_eco, loc_eco))
    p_down = sum(pmf_eco[x_eco < c_eco] * x_eco[x_eco < c_eco])
    p_eco = k_eco_m*(p_up + p_down)
    profit = p_biz + p_eco
    L = loss(c_biz, c_eco, k_biz_u, k_biz_m, k_eco_u, k_eco_m)
    net_profit = profit - L
    #print(L)
    NP[c_biz][c_eco] = net_profit

Максимальное значение средней общей прибыли составляет 70565 у.е. и достигается при $c^{(biz)}=19$ и $c^{(eco)}=81$ :

Python

np.max(NP)

Результат

70565.33152039198

Можно представить визуально уровень средней прибыли в зависимости от разных комбинаций мест:

Python

c_biz_max, c_eco_max = np.divmod(np.argmax(NP), 101)

im = plt.imshow(NP)
plt.xlim(65, 101)
plt.ylim(-1, 35)
plt.plot(c_eco_max, c_biz_max, 'ro', ms=4)
plt.vlines(c_eco_max, -1, c_biz_max, color='r', ls='--')
plt.hlines(c_biz_max, 0, c_eco_max, color='r', ls='--')
plt.text(70, c_biz_max - 2, r'$c^{*\mathrm {(biz)}}=19$')
plt.text(c_eco_max+1, 5, r'$c^{*\mathrm {(eco)}}=81$', ha='left')

plt.colorbar(im)
plt.ylabel('Number of seats in business class')
plt.xlabel('Number of seats in economy class')

plt.title('Dependence of the average total profit\n\
on the distribution of seats by class');

Значение не особо отличается от изначально выбранного авиакомпанией: 16 мест в бизнес классе и 84 места в экономе. Однако увеличение прибыли получается довольно значительным и составляет около 1.5%.

Python

100 * (NP[19, 81] - NP[16, 84]) / NP[16, 84]

Результат

1.4273417931471761

Если авиакомпания будет выбирать эти комбинации распределения кресел по классам случайно или основываясь на неверных или ошибочных суждениях, то потери потенциальной прибыли могут быть очень большими:

Python

lp = [100 * (NP[19, 81] - NP[i, 100 - i]) / NP[i, 100 - i] for i in range(10, 25)]

plt.plot(np.r_[10:25], lp, 'bo-', zorder=10)
plt.axhline(0, color='0.5')
plt.xlabel(r'$c^{biz}$ ($c^{eco} = N - c^{biz}$)')
plt.ylabel('%', rotation = 0, labelpad=10)
plt.title('Loss of average total profit');

В некоторой степени эта задача, как и предыдущая, все-таки оторвана от реальности. На самом деле авиакомпания не может просто взять и переделать внутреннюю конструкцию салона, ведь этим может заниматься только производитель. В данном случае речь может идти именно о покупке самолета, и лучше заранее знать оптимальное количество кресел в каждом подклассе перед заказом. Авиакомпания также может обратиться к производителю за переделкой салона, потому что текущая конфигурация ей невыгодна. При этом места в бизнес классе и эконом классе нельзя просто так суммировать, потому что в бизнес классе в одном ряду 4 кресла, а в экономе — 6. И даже такой нюанс довольно просто учесть в целевой функции.

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

Единственное, что сильно изменилось — это размерность. Может показаться, что дело пахнет керосином, ведь проклятие размерности никто не отменял. Для некоторых случаев это действительно так, но пока есть все предпосылки полагать, что значительная часть задач в отрасли гражданских авиаперевозок может быть решена обыкновенным брутфорсом (brute force, полный перебор). Та часть задач, которая не поддается брутфорсу, может же быть решена уже эвристическими методами. С интересным методом мы познакомимся далее, а пока рассмотрим способ извлечения еще большей выгоды от продажи авиабилетов.

Планирование продаж

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

Для увеличения доходов авиакомпании используют так называемое "классирование" пассажиров, которое обусловлено вполне естественным стремлением пассажиров к определенному уровню комфорта или экономии средств. Примерно так и появились всем знакомые классы обслуживания: первый, бизнес и эконом. Авиакомпании при этом пошли еще дальше, ведь внутри одного класса можно выделять разные уровни услуг. Например, можно взять бизнес класс и выделить в нем три разных уровня, назвав их:

МАКСИМУМ (FLEX);
ОПТИМУМ (CLASSIC);
БАЗОВЫЙ (BASE).

Данные уровни принято называть брендами. Бренды могут отличаться условиями предоставления услуг по перевозке: например, в бренде "МАКСИМУМ" можно вернуть билет без дополнительной платы, в бренде "ОПТИМУМ" с доплатой, а в бренде "БАЗОВЫЙ" возврат билета может быть вообще не разрешен. Условия предоставления услуги, как и её цена, прописываются в тарифе. Спрос на какой-нибудь отдельный бренд, допустим, на бренд "МАКСИМУМ" может постоянно меняться. Чтобы изменить цену билета соответствующим образом, необходимо изменить тариф. Это достаточно долгий процесс: например, в зарубежной системе ATPCO (Airline Tariff Publishing Company) изменения тарифа могут вступить в силу через час, а в отечественном ЦРТ (Центр Расписания и Тарифов) изменения вступят в силу только через сутки.

Чтобы менять стоимость билета на определенный бренд, авиакомпании используют классы бронирования, которые введены IATA (International Air Transport Association) как вспомогательный инструмент стандартизации условий путешествия пассажиров. Например, в бизнес классе существуют следующие классы бронирования:

J - Бизнес-класс, люкс (специальный тариф);
C - Бизнес-класс, полный тариф (нормальный /специальный тариф);
D - Бизнес-класс, льготный тариф (специальный тариф);
I - Бизнес-класс, льготный тариф (специальный тариф);
Z - Бизнес-класс, льготный тариф (специальный тариф).

Тарифы, привязанные к разным классам бронирования, могут иметь минимальные отличия в условиях оказания услуг: например, длительность остановки в пункте пересадки по тарифу класса бронирования J может быть не более 365 дней, а по тарифу класса бронирования C не более 360 дней. Главное отличие — это цена, и, чтобы менять её на некоторый бренд, достаточно менять продаваемый класс бронирования.

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

МАКСИМУМ (FLEX) - JFM, CFM, DFM, IFM, ZFM;
ОПТИМУМ (CLASSIC) - JCL, CCL, DCL, ICL, ZCL;
БАЗОВЫЙ (BASE) - JBS, CBS, DBS, IBS, ZBS.

Первая буква в коде тарифа означает класс бронирования. Тарифы с префиксом Z являются самыми дешевыми, а с J — самыми дорогими. При низком спросе можно продавать билеты по тарифам из класса бронирования Z, то есть ZFM, ZCL и ZBS. При увеличении спроса можно изменить класс бронирования продаваемых билетов, например с Z на D. После такого изменения купить билеты на соответствующий бренд будет можно только по более дорогим тарифам: DFM, DCL, DBS.

Пример выглядит странным, ведь закрытие одного класса бронирования и открытие другого означает, что стоимость билетов будет изменена сразу на все бренды. А если спрос меняется только на один бренд? В таком случае изменение класса бронирования приведет к оптимальному изменению цены по одному бренду и неоптимальному по другим. Увы, правила действительно так и работают. Максимум, что может придумать авиакомпания в таком случае — это выделить для каждого бренда хотя бы один уникальный класс бронирования, например:

МАКСИМУМ (FLEX) - JFM, IFM, ZFM;
ОПТИМУМ (CLASSIC) - CCL, ICL, ZCL;
БАЗОВЫЙ (BASE) - DBS, IBS, ZBS.

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

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

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

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

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

Отсюда возникает много вопросов к заявлениям как в научных статьях, так и в рекламной информации. Если какая-то компания (пусть даже очень авторитетная) заявляет, что её система автоматического ценообразования дает прирост прибыли в N %, то это вовсе не значит, что они действительно "выжимают" максимум. Методы ценообразования, как и данные о продажах, составляют коммерческую тайну, поэтому нет никакой возможности проверить адекватность используемых методов и реальный достигаемый максимум. Фактически верить авиакомпаниям приходится буквально на слово.

Некоторые методы все-таки есть: например, последние исследования (DPRank) показывают, что на основе исторических данных о ценах на билеты некоторой авиакомпании можно довольно точно восстановить информацию о спросе, то есть провести вполне адекватное сравнение методов ценообразования, используемых в авиакомпании с любыми другими.

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

Ограничимся теоретической моделью, поскольку я также не могу публиковать открытые данные по продажам. Рассмотрим 30-ти дневный период продаж билетов на некоторый бренд бизнес класса. Пусть величина спроса в некоторый день $\mathrm{day}$ будет задаваться распределением Пуассона, в котором параметр $\lambda$ также будет зависеть от количества дней до вылета:

$Q_{\mathrm{day}} \sim P(\lambda{(\mathrm{day})})$ $\lambda(\mathrm{day}) = \frac{a}{b + e^{c \cdot \mathrm{day} + d}} + k$

Функция $\lambda{(\mathrm{day})}$ представляет собой сигмоиду с дополнительными параметрами. Сигмоида является функцией, которая довольно хорошо описывает динамику спроса: сначала он долгое время остается неизменным, затем испытывает рост, а далее происходит его насыщение:

Python

def demand_by_day(n_day, ka, kb, kc, kd, ke):
    x = np.arange(0, n_day)
    return -ka / (kb + np.exp(-kc * x + kd)) + ke

q = demand_by_day(n_day=30, ka=3, kb=1.1, kc=0.35, kd=3, ke=4)
x = np.arange(0, 30)

plt.plot(x, q, 'bo', ms=2,
         label=r'$\lambda(\mathrm{day}) = \frac{3}{1.1 + e^{0.35 \cdot (\mathrm{day}) + 3}} + 4$')
plt.vlines(x, 0, q, colors='b', lw=1, alpha=0.5)
plt.legend()
plt.title(r'$\lambda(\mathrm{day})$')
plt.ylabel(r'$\lambda$ (persons)')
plt.xlabel('day');

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

Затем нам необходима кривая эластичности спроса по цене. Обычно в качестве кривой, описывающей эластичность, используются либо экспоненциальные, либо гиперболические функции. Как было сказано ранее, такие кривые показывают возможное количество покупателей по некоторой цене $\mathrm{price}$ . Поскольку мы вводим функцию, описывающую величину спроса в некоторый день, логично было бы описывать вероятность покупки билета (долю покупателей из ), то есть $P_{\mathrm{day}}(\mathrm{purchase} \; | \; \mathrm{price})$ . Тогда количество покупателей билета может быть описано с помощью биномиального распределения:

$B \sim \mathrm{Bin(Q, P_{\mathrm{day}}(\mathrm{purchase} \; | \; \mathrm{price}))}$

Для описания эластичности снова воспользуемся сигмоидой с дополнительными параметрами:

Python

pt_j = np.linspace(0, 10000, 500)
sd_j = -1 / (1 + np.exp(-0.001 * pt_j + 5)) + 1

plt.plot(pt_j, sd_j, 'b-', ms=3,
         label=r'$p = \frac{-1}{1 + e^{-0.001 \cdot (\mathrm{price}) + 5}} + 1$')
plt.legend()
plt.title('Price elasticity of demand model')
plt.ylabel(r'$P(\mathrm{purchase} \; | \; \mathrm{price})$')
plt.xlabel('price');

Теперь, когда нам доступна параметрическая модель, можно приступить к исследованию разных аспектов процесса продаж. Предположим, что авиакомпания продает билеты в бизнес класс вместимостью 40 кресел и использует для этого всего два класса бронирования: J и С, то есть всего две цены. Необходимо определить оптимальные значения цен: $\mathrm{price}_{J}$ и $\mathrm{price}_{C}$ , а также наилучший день $\mathrm{day^{*}}$ , в который происходит ее изменение:

$Z(\mathrm{price}_{J}^{*}, \mathrm{price}_{C}^{*}, \mathrm{day^{*}}) = \underset{\mathrm{price}_{J}, \mathrm{price}_{C}, \mathrm{day}}{\max} \left [ \mathrm{price}_{J} \cdot \! \sum_{\mathrm{day} = 1}^{\mathrm{day^{*}}} \!\! B_{\mathrm{day}} + \; \mathrm{price}_{C} \cdot \!\!\!\! \sum_{\mathrm{day} = \mathrm{day^{*}}}^{30} \!\!\! B_{\mathrm{day}} - \mathrm{Loss} \left ( N - \sum_{\mathrm{day} = 1}^{\mathrm{day^{*}}} \!\! B_{\mathrm{day}} - \sum_{\mathrm{day} = \mathrm{day^{*}}}^{30} \!\!\! B_{\mathrm{day}} \right ) \right ].$

Поскольку $B_{\mathrm{day}}$ — это количество проданных билетов в определенный день, а последнее место на рейсе должно быть продано в последний день продаж (день вылета), необходимо добавить ограничение:

$P(B_{0} > 0) > 0.95$

Данное ограничение необходимо, чтобы отсутствие билетов в последний день вылета наблюдалось не более чем в 5% случаев. Разумеется, хотелось бы вообще видеть 0% случаев, но к этому вопросу мы еще вернемся.

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

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

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

Смысл оптимизации методом кросс-энтропии заключается в повышении вероятности нахождения оптимального решения с помощью случайного перебора. Например, в рассматриваемой задаче вероятность нахождения оптимальных значений цен и дня их смены на первых порах перебора составит порядка 1:3000000000. Чтобы повысить вероятность нахождения оптимальных значений с помощью случайного перебора, достаточно правильным образом смещать матожидание и уменьшать дисперсию области, в которой выполняется перебор.

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

Python

# Функция, вычисляющая среднюю прибыль:
def e_profit(t, price_1, price_2, q):
    N = 40
    loss_es = 2500    # убытки от пустого кресла
    x = np.arange(0, 30)
    pd_1 = -1 / (1 + np.exp(-0.001 * price_1 + 5)) + 1
    pd_2 = -1 / (1 + np.exp(-0.001 * price_2 + 5)) + 1
    
    last_day = []
    Pr = []
    n_iter = 1000
    for i in range(n_iter):
        qd = poisson.rvs(mu=q)[::-1]
        n_bkd_1 = np.asarray(binom.rvs(n=qd[:t], p=pd_1))
        n_bkd_2 = np.asarray(binom.rvs(n=qd[t:], p=pd_2))
        n_bkd_1_sum, n_bkd_2_sum = n_bkd_1.sum(), n_bkd_2.sum()
        if n_bkd_1_sum <= N - n_bkd_2_sum:
            profit = price_1 * n_bkd_1_sum + price_2 * n_bkd_2_sum
        else:
            profit = price_1 * n_bkd_1_sum + price_2 * (N - n_bkd_1_sum)
        all_bkd = np.hstack((n_bkd_1, n_bkd_2))
        all_bkd_sum = all_bkd.sum()
        if all_bkd_sum < N:
            loss = (N - all_bkd_sum) * loss_es
            profit -= loss
        Pr.append(profit)
        last_day.append(N - np.sum(all_bkd[:29]))
    
    last_day = np.array(last_day)
    prob_up_0 = np.sum(last_day > 0) / n_iter
    return np.mean(Pr), prob_up_0

Python

# Функция, генерирующая случайные значения 
# переменных из заданных интервалов:
def vec_v(a_t, b_t, a_pr_1, b_pr_1, a_pr_2, b_pr_2):
    delta_p = 200
    if a_t != b_t:
        t = randint.rvs(a_t, b_t)
    else:
        t = a_t
    
    if a_pr_1 != b_pr_1:
        pr_1 = delta_p*(randint.rvs(a_pr_1, b_pr_1) // delta_p)
    else:
        pr_1 = a_pr_1
        
    if a_pr_2 != b_pr_2:
        pr_2 = delta_p*(randint.rvs(a_pr_2, b_pr_2) // delta_p)
    else:
        pr_2 = a_pr_2
    return t, pr_1, pr_2

Python

# Оптимизация методом кросс-энтропии:
'''
a_t, b_t = 1, 29
a_pr_1, b_pr_1 = 3000, 8000
a_pr_2, b_pr_2 = 3000, 8000

TIME, PRICE_1, PRICE_2 = [[a_t, b_t]], [[a_pr_1, b_pr_1]], [[a_pr_2, b_pr_2]]

for i in range(30):
    TIME.append([a_t, b_t])
    PRICE_1.append([a_pr_1, b_pr_1])
    PRICE_2.append([a_pr_2, b_pr_2])

    print(i)
    len_X = 0
    X, S = [], []

    while len_X < 100:
        x = vec_v(a_t, b_t, a_pr_1, b_pr_1, a_pr_2, b_pr_2)
        s, prob = e_profit(*x, q=q)
        if prob > 0.95 and s > 0:
            X.append(x)
            S.append(s)
            len_X += 1
            print('|', end='')
    X = np.array(X)
    S = np.array(S)
    idx = np.argsort(S)[::-1][:10]
    a_t, b_t = np.min(X[:, 0][idx]), np.max(X[:, 0][idx])
    a_pr_1, b_pr_1 = np.min(X[:, 1][idx]), np.max(X[:, 1][idx])
    a_pr_2, b_pr_2 = np.min(X[:, 2][idx]), np.max(X[:, 2][idx])

    TIME.append([a_t, b_t])
    PRICE_1.append([a_pr_1, b_pr_1])
    PRICE_2.append([a_pr_2, b_pr_2])
    
    print('\n', a_t, b_t, a_pr_1, b_pr_1, a_pr_2, b_pr_2, '\n'+'-'*20)
    if a_t == b_t and a_pr_1 == b_pr_1 and a_pr_2 == b_pr_2:
        break
''';

Благодаря данному методу можно довольно быстро найти оптимальное решение для представленной модели спроса:

$\mathrm{day^{*}}$ = 7;
$\mathrm{price}_{C}^{*}$ = 4400 у.е.
$\mathrm{price}_{J}^{*}$ = 5000 у.е.

Для извлечения максимальной выгоды необходимо 7 дней продавать билеты по тарифу из класса бронирования C то есть по цене 4400 у.е., а оставшиеся 23 дня билеты из J по 5000 у.е. Такая стратегия принесет в среднем 138600 у.е. прибыли:

Python

N = 40
t = 7
price_1 = 4400
price_2 = 5000
loss_es = 2500    # убытки от пустого кресла
pd_1 = -1 / (1 + np.exp(-0.001 * price_1 + 5)) + 1
pd_2 = -1 / (1 + np.exp(-0.001 * price_2 + 5)) + 1

last_day = []
Pr = []
Pr1 = []
Pr2 = []
LOSS = []
BKD_1 = []
BKD_2 = []
ALL = []
for i in range(50000):
    qd = poisson.rvs(mu=q)[::-1]
    n_bkd_1 = binom.rvs(n=qd[:t], p=pd_1)
    n_bkd_2 = binom.rvs(n=qd[t:], p=pd_2)
    if n_bkd_2.sum() <= N - n_bkd_1.sum():
        profit = price_1 * n_bkd_1.sum() + price_2 * n_bkd_2.sum()
        profit_1 = price_1 * n_bkd_1.sum()
        profit_2 = price_2 * n_bkd_2.sum()
    else:
        profit = price_1 * n_bkd_1.sum() + price_2 * (N - n_bkd_1.sum())
        profit_1 = price_1 * n_bkd_1.sum()
        profit_2 = price_2 * (N - n_bkd_1.sum())
    
    all_bkd = np.hstack((n_bkd_1, n_bkd_2))
    if all_bkd.sum() < N:
        loss = (N - all_bkd.sum()) * loss_es
        profit -= loss 
    Pr.append(profit)
    Pr1.append(profit_1)
    Pr2.append(profit_2)
    LOSS.append(loss)
    BKD_1.append(binom.rvs(n=qd[:t], p=pd_1).sum())
    BKD_2.append(binom.rvs(n=qd[t:], p=pd_2).sum())
    ALL.append(all_bkd.sum())
    last_day.append(40 - np.sum(all_bkd[:29]))

Python

np.mean(Pr)

Результат

138553.352

Мы также можем убедиться в том, что количество доступных билетов в последний день продаж наблюдается почти в 95% случаев:

Python

last_day = np.array(last_day)
b_leq = 100 * np.sum(last_day <= 0)/last_day.size
sns.histplot(last_day[last_day <= 0], discrete=True, label=f'{b_leq:.3}%')
b_lg = 100 * np.sum(last_day > 0)/last_day.size
sns.histplot(last_day[last_day > 0], discrete=True, label=f'{b_lg:.3}%')
plt.legend()
plt.text(-10, 2000, fr'$\sigma=${last_day.std():.3}')
plt.title('Number of tickets available on the last day of sales')
plt.xlabel('Number of tickets');

Теперь взглянем на общее количество продаваемых мест:

Python

ALL = np.array(ALL)
sns.histplot(ALL[ALL <= 40], discrete=True)
plt.title('Total number of tickets sold')
plt.axvline(ALL.mean(), color='red', label=f'{ALL.mean():.3}  tickets')
plt.legend()
plt.text(20, 2500, fr'$\sigma=${ALL[ALL <= 40].std():.3}')
plt.xlabel('Number of tickets');

Конечно, в идеале должно быть всегда 40, но давайте взглянем на распределение общего количества потенциальных покупателей:

Python

all_q = np.array([poisson.rvs(mu=q).sum() for i in range(10000)])
sns.histplot(all_q)
plt.title('Total number of potential buyers')
plt.xlabel('Number of potential buyers');

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

Можно ли с помощью стохастического программирования управлять заполняемостью мест? Для начала взглянем на количество продаваемых билетов в первые 7 дней продаж:

Python

sns.histplot(BKD_1, discrete=True)
plt.xticks(np.r_[0:18]);

Предположим, что количество проданных билетов, которое не входит в интервал [4, 8], считается сильным отклонением от приемлемой заполняемости и требует корректировки с помощью следующих цен:

если продано 0 билетов, то $\mathrm{price}_{J}$ = 4650 у.е.;
если продано 1 билетов, то $\mathrm{price}_{J}$ = 4750 у.е.;
если продано 2 билетов, то $\mathrm{price}_{J}$ = 4800 у.е.;
если продано 3 билетов, то $\mathrm{price}_{J}$ = 4900 у.е.;
если продано 9 билетов, то $\mathrm{price}_{J}$ = 5300 у.е.;
если продано 10 билетов, то $\mathrm{price}_{J}$ = 5350 у.е.;
если продано 11 билетов, то $\mathrm{price}_{J}$ = 5400 у.е.;
если продано 12 билетов, то $\mathrm{price}_{J}$ = 5450 у.е.;
если продано 13 билетов, то $\mathrm{price}_{J}$ = 5550 у.е.;
если продано 14 билетов, то $\mathrm{price}_{J}$ = 5600 у.е.;
если продано 15 билетов, то $\mathrm{price}_{J}$ = 5700 у.е..

Вычислить корректирующие цены для данного примера можно с помощью простого перебора:

Python

'''t = 7
B_1 = np.hstack((np.r_[0:4], np.r_[9:18]))
cor_prices = []
for b_1 in B_1:
    N = 40 - b_1
    price_1 = 4400
    price_2 = 5000
    loss_es = 2500    # убытки от пустого кресла
    x = np.arange(0, 30)
    pd_1 = -1 / (1 + np.exp(-0.001 * price_1 + 5)) + 1
    pd_2 = -1 / (1 + np.exp(-0.001 * price_2 + 5)) + 1
    BKD_1 = []
    BKD_2 = []
    last_day = []
    Pr = []
    n_iter = 4000
    PROFIT = []
    m_profit = 0
    pr_max = 0
    pr_x = np.arange(3900, 8000, 50)
    for price_2 in pr_x:
        #print(price_2)
        last_day = []
        Pr = []
        for i in range(n_iter):
            qd = poisson.rvs(mu=q)[::-1]
            n_bkd_1 = b_1
            pd_2 = -1 / (1 + np.exp(-0.001 * price_2 + 5)) + 1
            n_bkd_2 = np.asarray(binom.rvs(n=qd, p=pd_2))
            last_day.append(N - n_bkd_2[:-1].sum())
            n_bkd_2 = n_bkd_2.sum()
            BKD_2.append(n_bkd_2)
            if n_bkd_2 <= N:
                profit = price_1 * n_bkd_1 + price_2 * n_bkd_2
            else:
                profit = price_1 * n_bkd_1 + price_2 * N

            if n_bkd_2 < N:
                loss = (N - n_bkd_2) * loss_es
                profit -= loss 
            Pr.append(profit)

        last_day = np.array(last_day)
        prob = np.sum(last_day > 0) / n_iter
        PROFIT.append(np.mean(Pr))
        if prob >= 0.95:
            if m_profit < np.mean(Pr):
                m_profit = np.mean(Pr)
                pr_max = price_2
                break
    print(b_1, pr_max, np.sum(last_day > 0) / n_iter)
    cor_prices.append(pr_max)
cor_prices'''

# результат:
cor_prices = [4650, 4700, 4800, 4900, 5300, 5350, 5400, 5500, 5500, 5600, 5700, 5750, 5850]

Python

N = 40
t = 7
price_1 = 4400
pd_1 = -1 / (1 + np.exp(-0.001 * price_1 + 5)) + 1
loss_es = 2500    # убытки от пустого кресла


last_day = []
Pr_star = []
Pr1 = []
Pr2 = []
LOSS = []
BKD_1 = []
BKD_2 = []
ALL = []
Price_2 = []
for i in range(50000):
    qd = poisson.rvs(mu=q)[::-1]
    n_bkd_1 = binom.rvs(n=qd[:t], p=pd_1)
    n_bkd_1_sum = n_bkd_1.sum()
    if n_bkd_1_sum < 4 or n_bkd_1_sum > 8:
        if n_bkd_1_sum == 0:
            price_2 = 4650
        if n_bkd_1_sum == 1:
            price_2 = 4750
        if n_bkd_1_sum == 2:
            price_2 = 4800
        if n_bkd_1_sum == 3:
            price_2 = 4900
        if n_bkd_1_sum == 9:
            price_2 = 5300
        if n_bkd_1_sum == 10:
            price_2 = 5350
        if n_bkd_1_sum == 11:
            price_2 = 5400
        if n_bkd_1_sum == 12:
            price_2 = 5450
        if n_bkd_1_sum == 13:
            price_2 = 5550
        if n_bkd_1_sum == 14:
            price_2 = 5600
        if n_bkd_1_sum == 15:
            price_2 = 5700
        if n_bkd_1_sum == 16:
            price_2 = 5750
        if n_bkd_1_sum == 17:
            price_2 = 5850
    else:
        price_2 = 5000
    Price_2.append(price_2)
    pd_2 = -1 / (1 + np.exp(-0.001 * price_2 + 5)) + 1
    n_bkd_2 = binom.rvs(n=qd[t:], p=pd_2)
    all_bkd = np.hstack((n_bkd_1, n_bkd_2))
    if n_bkd_2.sum() <= N - n_bkd_1.sum():
        profit = price_1 * n_bkd_1.sum() + price_2 * n_bkd_2.sum()
        profit_1 = price_1 * n_bkd_1.sum()
        profit_2 = price_2 * n_bkd_2.sum()
    else:
        profit = price_1 * n_bkd_1.sum() + price_2 * (N - n_bkd_1.sum())
        profit_1 = price_1 * n_bkd_1.sum()
        profit_2 = price_2 * (N - n_bkd_1.sum())
    if all_bkd.sum() < N:
        loss = (N - all_bkd.sum()) * loss_es
        profit -= loss 
    Pr_star.append(profit)
    Pr1.append(profit_1)
    Pr2.append(profit_2)
    LOSS.append(loss)
    BKD_1.append(binom.rvs(n=qd[:t], p=pd_1).sum())
    BKD_2.append(binom.rvs(n=qd[t:], p=pd_2).sum())
    ALL.append(n_bkd_1_sum + n_bkd_2.sum())
    last_day.append(N - np.sum(all_bkd[:29]))

Применение корректировок не привело к снижению средней прибыли, и она по-прежнему составляет 138600 у.е.

Python

np.mean(Pr)

Результат

138553.352

В то же время, количество доступных билетов в последний день продаж наблюдается на 2% чаще, чем без корректировок, а дисперсия этого количества уменьшилась более чем на 7%:

Python

last_day = np.array(last_day)
b_leq = 100 * np.sum(last_day <= 0)/last_day.size
sns.histplot(last_day[last_day <= 0], discrete=True, label=f'{b_leq:.3}%')
b_lg = 100 * np.sum(last_day > 0)/last_day.size
sns.histplot(last_day[last_day > 0], discrete=True, label=f'{b_lg:.3}%')
plt.legend()
plt.text(-10, 2000, fr'$\sigma=${last_day.std():.3}')
plt.title('Number of tickets available on the last day of sales')
plt.xlabel('Number of tickets');

В данном случае мы выполняли корректировку для удовлетворения изначально установленному ограничению:

$P(B_{0} > 0) > 0.95$

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

Python

price_1 = 4400
price_2 = 5000

pd_1 = -1 / (1 + np.exp(-0.001 * price_1 + 5)) + 1
pd_2 = -1 / (1 + np.exp(-0.001 * price_2 + 5)) + 1
good_occupancy = []
bad_occupancy = []
for i in range(3000):
    qd = poisson.rvs(mu=q)[::-1]
    n_bkd_1 = binom.rvs(n=qd[:t], p=pd_1)
    n_bkd_2 = binom.rvs(n=qd[t:], p=pd_2)
    all_bkd = np.cumsum(np.hstack((n_bkd_1, n_bkd_2)))
    if 39 <= all_bkd[-1] <= 40:
        good_occupancy.append(all_bkd)
    else:
        bad_occupancy.append(all_bkd)
        
good_occupancy_ser = pd.Series(np.hstack(good_occupancy),
                           index=list(range(30))*len(good_occupancy))

Python

for occ in good_occupancy:
    plt.plot(occ, 'bo', ms=1, alpha=0.2, zorder=10)
sns.lineplot(good_occupancy_ser, errorbar=('pi', 80))
plt.title('optimal filling process')
plt.ylabel('Number of tickets')
plt.xlabel('day');

На данном графике показан 80% процентильный интервал. Если в какой-то день количество проданных мест выходит за пределы данного интервала, то можно подобрать цену, которая с наибольшей вероятностью вернет процесс заполняемости в данный интервал через какое-то фиксированное количество дней. Именно это и будет являться задачей динамического ценообразования.

Благодаря планированию продаж возможно ответить на многие важные вопросы. Например, если кривая эластичности с течением времени испытывает сдвиги вправо или влево, то можно обнаружить, что либо $\mathrm{price}_{C}$ станет больше, чем $\mathrm{price}_{J}$ , либо $\mathrm{day^{*}}$ станет равен 30 дням до вылета. В первом случае, когда спрос большой, оптимальнее будет начинать продажи за меньшее количество дней до вылета, причем начинать со дня $\mathrm{day^{*}}$ . В случае, когда спрос невелик и $\mathrm{day^{*}}$ =30, то день начала продаж необходимо отодвигать на более ранние сроки. В зависимости от ограничений такое поведение может наблюдаться и при других изменениях спроса.

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

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

Python

def e_profit(J, C, price_1, price_2, q):
    N = 40
    loss_es = 2500    # убытки от пустого кресла

    pd_1 = -1 / (1 + np.exp(-0.001 * price_1 + 5)) + 1
    pd_2 = -1 / (1 + np.exp(-0.001 * price_2 + 5)) + 1

    last_day = []
    Pr = []
    n_iter = 30000

    for i in range(n_iter):
        qd = poisson.rvs(mu=q)[::-1]
        bkd_1 = np.asarray(binom.rvs(n=qd, p=pd_2))
        diff_1 = bkd_1.cumsum()-C

        t1 = diff_1[diff_1 < 0].size

        if t1 < 30:
            bkd_1 = bkd_1[:t1 + 1].sum()
            if bkd_1 >= C:
                bkd_1 = C
            qd[t1] -= diff_1[t1]

            bkd_2 = np.asarray(binom.rvs(n=qd[t1:], p=pd_1))
            diff_2 = bkd_2.cumsum()-J

            t2 = diff_2[diff_2 < 0].size

            bkd_2 = bkd_2.sum()
            if bkd_2 >= J:
                bkd_2 = J

        else:
            bkd_1 = bkd_1[:t1 + 1].sum()
            if bkd_1 >= C:
                bkd_1 = C
            t2 = 30
            bkd_2 = 0

        profit = price_2 * bkd_1 + price_1 * bkd_2 - loss_es * (N - bkd_1 - bkd_2)

        Pr.append(profit)
        last_day.append(t1+t2)
    last_day = np.array(last_day)
    prob = last_day[last_day == 30].size / n_iter
    return np.mean(Pr), prob

def vec_v(a_j, b_j, a_pr_1, b_pr_1, a_pr_2, b_pr_2):
    if a_j == b_j:
        quota_j = a_j
        quota_c = 40 - quota_j
    else:
        quota_j = randint.rvs(low=a_j, high=b_j)
        quota_c = 40 - quota_j
    
    delta_p = 200
    if a_pr_1 != b_pr_1:
        pr_1 = delta_p*(randint.rvs(a_pr_1, b_pr_1) // delta_p)
    else:
        pr_1 = a_pr_1
        
    if a_pr_2 != b_pr_2:
        pr_2 = delta_p*(randint.rvs(a_pr_2, b_pr_2) // delta_p)
    else:
        pr_2 = a_pr_2
    return quota_j, quota_c, pr_1, pr_2

'''
a_j, b_j = 20, 39
a_pr_1, b_pr_1 = 4000, 7000
a_pr_2, b_pr_2 = 4000, 7000

for i in range(10):
    print(i)
    len_X = 0
    X, S = [], []

    while len_X < 100:
        x = vec_v(a_j, b_j, a_pr_1, b_pr_1, a_pr_2, b_pr_2)
        s, prob = e_profit(*x, q=q)
        if prob > 0.95 and x[2] < x[3]:
            X.append(x)
            S.append(s)
            len_X += 1
            print('|', end='')
    X = np.array(X)
    S = np.array(S)
    idx = np.argsort(S)[::-1][:20]
    a_j, b_j = np.min(X[:, 0][idx]), np.max(X[:, 0][idx])
    a_pr_1, b_pr_1 = np.min(X[:, 2][idx]), np.max(X[:, 2][idx])
    a_pr_2, b_pr_2 = np.min(X[:, 3][idx]), np.max(X[:, 3][idx])
    print('\n', a_j, b_j, a_pr_1, b_pr_1, a_pr_2, b_pr_2, '\n'+'-'*20)
'''
# оптимальное решение: 29, 11, 5000, 5400

Для данной модели спроса оптимальная квота мест в классе бронирования J составит 29 мест по цене 5000 у.е., а для C — 11 мест по 5400 у.е. Выглядит нелогично, но это объясняется тем, что модель ищет наилучший способ удовлетворения ограничению. Если использование квот является принципиально необходимым, то можно указать дополнительное ограничение $\mathrm{price}_{C} < \mathrm{price}_{J}$ . Средняя прибыль в такой стратегии значительно меньше и составит всего 123600 у.е. Использование квот приводит к небольшому уменьшению матожидания и небольшому увеличению отклонения общего количества продаваемых билетов, что еще раз говорит о неочевидном поведении арифметических действий над случайными величинами.

Python

N = 40
loss_es = 2500    # убытки от пустого кресла
price_1, price_2 = 5000, 5400
C, J = 11, 29
pd_1 = -1 / (1 + np.exp(-0.001 * price_1 + 5)) + 1
pd_2 = -1 / (1 + np.exp(-0.001 * price_2 + 5)) + 1

last_day = []
Pr = []
BKD_2 = []
n_iter = 100000

for i in range(n_iter):
    qd = poisson.rvs(mu=q)[::-1]
    bkd_1 = np.asarray(binom.rvs(n=qd, p=pd_2))
    diff_1 = bkd_1.cumsum()-C

    t1 = diff_1[diff_1 < 0].size

    if t1 < 30:
        bkd_1 = bkd_1[:t1 + 1].sum()
        if bkd_1 >= C:
            bkd_1 = C
        qd[t1] -= diff_1[t1]

        bkd_2 = np.asarray(binom.rvs(n=qd[t1:], p=pd_1))
        diff_2 = bkd_2.cumsum()-J

        t2 = diff_2[diff_2 < 0].size

        bkd_2 = bkd_2.sum()
        if bkd_2 >= J:
            bkd_2 = J

    else:
        bkd_1 = bkd_1[:t1 + 1].sum()
        if bkd_1 >= C:
            bkd_1 = C
        t2 = 30
        bkd_2 = 0

    profit = price_2 * bkd_1 + price_1 * bkd_2 - loss_es * (N - bkd_1 - bkd_2)

    Pr.append(profit)
    last_day.append(t1+t2)
    BKD_2.append(bkd_2)

Python

sns.histplot(np.array(BKD_2) + 11, discrete=True)
plt.title('Total number of tickets sold')
plt.axvline(np.mean(np.array(BKD_2) + 11), color='red',
            label=f'{np.mean(np.array(BKD_2) + 11):.3}  tickets')
plt.legend()
plt.text(15, 2500, fr'$\sigma=${np.std(np.array(BKD_2) + 11):.3}')
plt.xlabel('Number of tickets');

Разумеется, небольшие отклонения от оптимальных значений могут приводить к сильным отклонениям средней прибыли. Например, для изначальной самой прибыльной стратегии увеличение $\mathrm{price}_{J}$ всего на 1% обернется уменьшением средней прибыли на 2.5%. Если уменьшить $\mathrm{price}_{J}$ , то это приведет к увеличению прибыли, но и установленное ограничение будет нарушено:

Python

'''
res_pr = []
for price_j in range(4850, 5151, 20):
    res_pr.append(e_profit(7, 4400, price_j, q=q))

res_pr = np.array(res_pr)
'''
res_pr = np.array([[1.46966595e+05, 8.89166667e-01],
       [1.46322563e+05, 8.96666667e-01],
       [1.44654792e+05, 9.06066667e-01],
       [1.43902872e+05, 9.11900000e-01],
       [1.42776907e+05, 9.20000000e-01],
       [1.41942047e+05, 9.24700000e-01],
       [1.40439868e+05, 9.30933333e-01],
       [1.39520810e+05, 9.39933333e-01],
       [1.37781890e+05, 9.42633333e-01],
       [1.36935494e+05, 9.47333333e-01],
       [1.35520077e+05, 9.53566667e-01],
       [1.33495468e+05, 9.60233333e-01],
       [1.32832887e+05, 9.61800000e-01],
       [1.30971382e+05, 9.65700000e-01],
       [1.30005971e+05, 9.69766667e-01],
       [1.28637470e+05, 9.72700000e-01]])

percent = 100 * (res_pr[:, 0] - 138600) / 138600
plt.plot(percent, 'g-')
plt.plot(np.r_[0:16][res_pr[:, 1] < 0.94],
         percent[res_pr[:, 1] < 0.94], 'rX')
plt.plot(np.r_[0:16][res_pr[:, 1] > 0.94],
         percent[res_pr[:, 1] > 0.94], 'bo')
plt.axhline(0, color='0.4')
plt.xticks(np.r_[0:16], range(4850, 5151, 20), rotation=45);

Еще следует сделать небольшое отступление про метод оптимизации с помощью кросс-энтропии. Это действительно очень мощный инструмент, который позволяет подобраться к, казалось бы, неприступным задачам. Теоретически EM (Enterprise management) и RM как дисциплины решают только одну задачу, но на практике оказывается, что её необходимо разбивать на подзадачи. В определенной степени умение выделять подзадачи считается искусством.

Например, задача управления авиакомпанией разбивается на подзадачи:

поиск оптимальной маршрутной сети;
назначении флота и экипажей;
составление расписания и т.д.

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

Оптимизация с помощью кросс-энтропии позволяет объединять задачи. В результате объединения получаются задачи, к которым поначалу даже непонятно как подступиться, например, комбинаторно-стохастические задачи очень характерны для логистических цепочек. Есть и множество других сложных задач, например, о справедливом разрезании торта.

Все это совершенно точно не является готовым рецептом, это скорее концепция. И ее адаптация к определенной задаче может быть по-своему трудным процессом. Например, используемый для демонстрации вариант алгоритма является очень "неосторожным", ведь изначальное ограничение по наличию свободных мест в 95% случаев не было соблюдено в точности. В полученном решении оно гарантировало 94%. Это значит, что на каждой итерации необходимо выполнять дополнительные статистические тесты. Также алгоритм является очень медленным: чтобы добиться гораздо большего прироста в скорости, необходимо пользоваться другими распределениями. Используемый алгоритм также можно назвать "расточительным", ведь множество точек генерируется заново на каждой итерации вместо того, чтобы отсеивать неподходящие и добавлять новые для анализа.

Адаптировать данный метод оптимизации к задачам можно итеративно, а его важным преимуществом является то, что он позволяет исследовать пространство возможных решений. Нужно понимать, что такие алгоритмы, как симуляция отжига или рой частиц, являются подходящими только для определенных классов задач. На практике появляются задачи (модели), которые без преувеличения можно назвать особенными — именно для таких задач оптимизация с помощью кросс-энтропии является наилучшим методом.

Проблема отслеживания динамики спроса

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

Создание модели спроса — это стандартная задача аппроксимации, для решения которой можно использовать самые разные способы: регрессия, случайный лес, нейронные сети или бустинг. Есть достаточно оснований полагать, что для анализа спроса регрессия на основе гаусовских процессов подходят лучше всего. Однако не все так просто, как хотелось бы, ведь спрос зависит не только от цены, но еще от других факторов: наличия альтернативных предложений и цен на них, погоды, дня недели, сезона и т.д. Чем больше факторов учитывается — тем лучше. В RM речь идет о "здесь и сейчас", а значит актуальных данных мало, к тому же цены на самом деле ограничивают и цензурируют наблюдения. А еще ведь бронирования бывают групповые! Бывают возвраты и обмены билетов, так что никакой ярко выраженной "нормальности" в данных нет, и для их предварительной обработки необходимо прикладывать немало усилий.

Модель спроса является это типичным примером модели стационарной системы, которая пригодна для того, чтобы понять, как что-то происходит в "среднем по больнице". Динамическое ценообразование на основе таких моделей не может быть полностью автоматизировано, поскольку оно не учитывает факторы, которые невозможно заложить в модель: техногенные и природные катастрофы, политические действия разных стран, резкие скачки в экономике и т.д. Такие системы могут отслеживать резкие отклонения спроса и просигнализировать о необходимости дополнительного анализа ситуации, но реагировать на резкие отклонения самостоятельно они не способны.

В одной из предыдущих статей мы писали об отслеживании динамики с помощью критерия Вальда. В рассматриваемой модели речь шла о всей емкости довольно большого самолета, и этого небольшого количества наблюдений было достаточно, чтобы судить об изменениях спроса с приемлемыми уровнями ошибок. В данном же случае мы рассматриваем отдельные классы обслуживания и классы бронирования. Это значит, что объемы спроса могут быть недостаточны для критерия Вальда — нужно что-то более чувствительное.

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

Например, возможные переходы количества купленных билетов в 6-ом и 7-ом дне продаж для представленной модели спроса будут выглядеть следующим образом:

Python

price_1 = 4400
price_2 = 5000

pd_1 = -1 / (1 + np.exp(-0.001 * price_1 + 5)) + 1
pd_2 = -1 / (1 + np.exp(-0.001 * price_2 + 5)) + 1
transitions = []
n_iter = 10000
for i in range(n_iter):
    qd = poisson.rvs(mu=q)[::-1]
    n_bkd_1 = binom.rvs(n=qd[[5, 6]], p=pd_1)
    n_bkd_2 = binom.rvs(n=qd[7], p=pd_2)
    transitions.append(np.cumsum(np.hstack((n_bkd_1, n_bkd_2))))
    
transitions = np.array(transitions)

transitions_6_7, count_6_7 = np.unique(transitions[:, :2],
                                       axis=0, return_counts= True)
transitions_7_8, count_7_8 = np.unique(transitions[:, 1:],
                                       axis=0, return_counts= True)
idx_sort_6_7 = np.argsort(count_6_7)
idx_sort_7_8 = np.argsort(count_7_8)

transitions_6_7, count_6_7 = transitions_6_7[idx_sort_6_7], count_6_7[idx_sort_6_7]
transitions_7_8, count_7_8 = transitions_7_8[idx_sort_7_8], count_7_8[idx_sort_7_8]

probs_6_7 = count_6_7 / n_iter
probs_7_8 = count_7_8 / n_iter

DG = nx.DiGraph()
DG.add_edges_from(transitions_6_7)
plt.figure(figsize=(5, 5))
pos = nx.circular_layout(DG)
#nx.draw_networkx(DG, pos=pos)
nx.draw_networkx_nodes(DG, pos=pos, node_color='#ffffe4',
                 edgecolors='0.4')
nx.draw_networkx_labels(DG, pos=pos)
edges_red = nx.draw_networkx_edges(DG, pos=pos, 
                   edgelist=transitions_6_7[probs_6_7.cumsum() < 0.05],
                   edge_color='#cb416b')
edges_red[0].set_label('p = 0.05')
edges_blue = nx.draw_networkx_edges(DG, pos=pos, 
                   edgelist=transitions_6_7[probs_6_7.cumsum() > 0.05],
                   edge_color='#601ef9')
edges_blue[0].set_label('p = 0.95')
plt.legend()
plt.title('Possible transitions in sales volumes\non the 6th and 7th days of sales');

Красные дуги означают переходы, которые для рассматриваемой модели спроса будут происходить лишь в 5% случаев. Если наблюдается несколько таких событий или событий, которых даже нет в цепи, то это становится сигналом к перерасчету параметров модели, а также оптимальных цен. К сожалению, "на лету" возможен только грубый перерасчет параметров, но главное, что такая возможность нам вообще доступна. Ведь, получив какое-то аномальное событие, можно найти минимальные отклонения параметров, которые могли бы привести к такому событию — именно такой перерасчет становится отдельной оптимизационной задачей.

Заключение

Представленные примеры показывают, что потери прибыли могут исчисляться целыми десятками процентов. В лучшем случае авиакомпания может обладать одной из предлагаемых на рынке RMS, при этом нам никогда и ничего неизвестно, что находится "под капотом" этих RMS. С одной стороны, это естественно, ведь для разработчиков RMS важно, чтобы их наработки охранялись коммерческой тайной. С другой стороны, в тайне можно хранить не только передовые наработки, но и устаревшие, неоптимальные методы. Также важно учитывать, что закрытость может вести к стагнации и последующей деградации любой отрасли.

Тема RM очень обширна, и в данной статье пока удалось охватить лишь несколько частей данной дисциплины — оптимизацию распределения ресурсов и планирование продаж. Тем не менее, по полученным результатам уже можно сделать вывод, что поиск оптимальных решений в условиях неопределенности возможен. Хочется верить, что данный материал рано или поздно поможет авиакомпаниям подходить более тщательно к RM или выбору специализированного ПО и даст мотивацию к изучению и правильному продвижению методов RM в других бизнесах, отраслях и сферах. Я же в одной из будущих статей хочу развить тему дальше и рассмотреть нестинг в части управления доходами.

Логистика. Часть 5. Управление доходами, или первый шаг к нестингу