根据 Python mystic 中的条件最大化总和
Maximize sum subject to conditions in Python mystic
我正在尝试在 Dynamic room pricing model for hotel revenue management systems 上构建白皮书的实现。以防这个 link 将来死亡，我将粘贴到此处的相关部分：
到目前为止，我目前的实现非常糟糕，因为我真的没有完全理解如何求解非线性最大化方程。
# magical lookup table that returns demand based on those inputs
# this will eventually be a db lookup against past years rental activity and not hardcoded to a specific value
def demand(dateFirstNight, duration):
return 1
# magical function that fetches the price we have allocated for a room on this date to existing customers
# this should be a db lookup against previous stays, and not hardcoded to a specific value
def getPrice(date):
return 75
# Typical room base price
# Defined as: Nominal price of the hotel (usually the average historical price)
nominalPrice = 89
# from the white paper, but perhaps needs to be adjusted in the future using the methods they explain
priceElasticity = 2
# this is an adjustable constant it depends how far forward we want to look into the future when optimizing the prices
# likely this will effect how long this will take to run, so it will be a balancing game with regards to accuracy vs runtime
numberOfDays = 30
def roomsAlocated(dateFirstNight, duration)
roomPriceSum = 0.0
for date in range(dateFirstNight, dateFirstNight+duration-1):
roomPriceSum += getPrice(date)
return demand(dateFirstNight, duration) * (roomPriceSum/(nominalPrice*duration))**priceElasticity
def roomsReserved(date):
# find all stays that contain this date, this
def maximizeRevenue(dateFirstNight):
# we are inverting the price sum which is to be maximized because mystic only does minimization
# and when you minimize the inverse you are maximizing!
return (sum([getPrice(date)*roomsReserved(date) for date in range(dateFirstNight, dateFirstNight+numberOfDays)]))**-1
def constraint(x): # Ol - totalNumberOfRoomsInHotel <= 0
return roomsReserved(date) - totalNumberOfRoomsInHotel
from mystic.penalty import quadratic_inequality
@quadratic_inequality(constraint, k=1e4)
def penalty(x):
return 0.0
from mystic.solvers import diffev2
from mystic.monitors import VerboseMonitor
mon = VerboseMonitor(10)
bounds = [0,1e4]*numberOfDays
result = diffev2(maximizeRevenue, x0=bounds, penalty=penalty, npop=10, gtol=200, disp=False, full_output=True, itermon=mon, maxiter=M*N*100)
任何熟悉使用 mystic 的人都可以给我一些关于如何实现它的建议吗？
虽然您要求使用库 mystic，但在开始非线性优化时，您可能不需要这种细粒度的控制。模块 scipy 应该足够了。以下是一个或多或少的完整解决方案，更正了我认为原始白皮书中关于定价范围的错字：
import numpy as np
from scipy.optimize import minimize
P_nom = 89
P_max = None
price_elasticity = 2
number_of_days = 7
demand = lambda a, L: 1./L
total_rooms = [5]*number_of_days
def objective(P, *args):
return -np.dot(P, O(P, *args))
def worst_leftover(P, C, *args):
return min(np.subtract(C, O(P, *args)))
def X(P, a, L, d, e, P_nom):
return d(a, L)*(sum(P[a:a+L])/P_nom/L)**e
def d(a, L):
return 1.
def O_l(P, l, l_max, *args):
return sum([X(P, a, L, *args)
for a in xrange(0, l)
for L in xrange(l-a+1, l_max+1)])
def O(P, *args):
return [O_l(P, l, *args) for l in xrange(len(P))]
result = minimize(
objective,
[P_nom]*number_of_days,
args=(number_of_days-1, demand, price_elasticity, P_nom),
method='SLSQP',
bounds=[(0, P_max)]*number_of_days,
constraints={
'type': 'ineq',
'fun': worst_leftover,
'args': (total_rooms, number_of_days-1, demand, price_elasticity, P_nom)
},
tol=1e-1,
options={'maxiter': 10**3}
)
print result.x
有几点值得一提：
objective 函数添加了一个减号，用于 scipy 的 minimize() 例程，与原始白皮书中引用的最大化形成对比。这将导致 result.fun 为负数而不是表示总收入。
公式好像对参数有点敏感。最小化是正确的（至少，当它说它正确执行时它是正确的——检查 result.success），但如果输入与现实相差太远，那么您可能会发现价格比预期高得多。您可能还想使用比我在以下示例中使用的天数更多的天数。您的白皮书似乎会产生周期性影响。
我不太喜欢白皮书的命名方案，因为它与可读代码有关。我改变了一些东西，但有些东西确实很糟糕，应该被替换，比如小写的 l 很容易与数字 1 混淆。
我确实设置了界限，使价格为正而不是负。凭借您的领域专业知识，您应该验证这是正确的决定。
您可能更喜欢比我指定的更严格的公差。这在某种程度上取决于您希望运行时是什么。随意使用 tol 参数。此外，对于更严格的公差，您可能会发现 options 参数中的 'maxiter' 必须增加才能使 minimize() 正确收敛。
我很确定 total_rooms 应该是酒店中尚未预订的房间数量，因为白皮书将其索引为字母 l 而不是常量你有你的原始代码。我将其设置为用于测试目的的常量列表。
方法必须是 'SLSQP' 来处理价格范围和房间数量范围。注意不要改这个。
O_l() 的计算效率非常低下。如果运行时是个问题，我要采取的第一步是弄清楚如何 cache/memoize 调用 X()。所有这些实际上只是第一步，概念验证。它应该没有错误且正确，但它几乎是直接从白皮书中提取的，可以进行一些重构。
任何人，玩得开心，如有任何其他问题，请随时 comment/PM/etc。
抱歉，我来晚了，但我认为接受的答案并没有解决完整的问题，而且还错误地解决了问题。请注意，在局部最小化中，求解接近名义价格的问题并不能给出最佳解决方案。
我们先建一个hotel class:
"""
This file is 'hotel.py'
"""
import math
class hotel(object):
def __init__(self, rooms, price_ave, price_elastic):
self.rooms = rooms
self.price_ave = price_ave
self.price_elastic = price_elastic
def profit(self, P):
# assert len(P) == len(self.rooms)
return sum(j * self._reserved(P, i) for i,j in enumerate(P))
def remaining(self, P): # >= 0
C = self.rooms
# assert len(P) == C
return [C[i] - self._reserved(P, i) for i,j in enumerate(P)]
def _reserved(self, P, day):
max_days = len(self.rooms)
As = range(0, day)
return sum(self._allocated(P, a, L) for a in As
for L in range(day-a+1, max_days+1))
def _allocated(self, P, a, L):
P_nom = self.price_ave
e = self.price_elastic
return math.ceil(self._demand(a, L)*(sum(P[a:a+L])/(P_nom*L))**e)
def _demand(self, a,L): #XXX: fictional non-constant demand function
return abs(1-a)/L + 2*(a**2)/L**2
这是使用 mystic 解决问题的一种方法：
"""
This file is 'local.py'
"""
n_days = 7
n_rooms = 50
P_nom = 85
P_bounds = 0,None
P_elastic = 2
import hotel
h = hotel.hotel([n_rooms]*n_days, P_nom, P_elastic)
def objective(price, hotel):
return -hotel.profit(price)
def constraint(price, hotel): # <= 0
return -min(hotel.remaining(price))
bounds = [P_bounds]*n_days
guess = [P_nom]*n_days
import mystic as my
@my.penalty.quadratic_inequality(constraint, kwds=dict(hotel=h))
def penalty(x):
return 0.0
# using a local optimizer, starting from the nominal price
solver = my.solvers.fmin
mon = my.monitors.VerboseMonitor(100)
kwds = dict(disp=True, full_output=True, itermon=mon,
args=(h,), xtol=1e-8, ftol=1e-8, maxfun=10000, maxiter=2000)
result = solver(objective, guess, bounds=bounds, penalty=penalty, **kwds)
print([round(i,2) for i in result[0]])
结果：
>$ python local.py
Generation 0 has Chi-Squared: -4930.000000
Generation 100 has Chi-Squared: -22353.444547
Generation 200 has Chi-Squared: -22410.402420
Generation 300 has Chi-Squared: -22411.780268
Generation 400 has Chi-Squared: -22413.908944
Generation 500 has Chi-Squared: -22477.869093
Generation 600 has Chi-Squared: -22480.144244
Generation 700 has Chi-Squared: -22480.280379
Generation 800 has Chi-Squared: -22485.563188
Generation 900 has Chi-Squared: -22485.564265
Generation 1000 has Chi-Squared: -22489.341602
Generation 1100 has Chi-Squared: -22489.345912
Generation 1200 has Chi-Squared: -22489.351219
Generation 1300 has Chi-Squared: -22491.994305
Generation 1400 has Chi-Squared: -22491.994518
Generation 1500 has Chi-Squared: -22492.061127
Generation 1600 has Chi-Squared: -22492.573672
Generation 1700 has Chi-Squared: -22492.573690
Generation 1800 has Chi-Squared: -22492.622064
Generation 1900 has Chi-Squared: -22492.622230
Optimization terminated successfully.
Current function value: -22492.622230
Iterations: 1926
Function evaluations: 3346
STOP("CandidateRelativeTolerance with {'xtol': 1e-08, 'ftol': 1e-08}")
[1.15, 20.42, 20.7, 248.1, 220.56, 41.4, 160.09]
再次使用全局优化器：
"""
This file is 'global.py'
"""
n_days = 7
n_rooms = 50
P_nom = 85
P_bounds = 0,None
P_elastic = 2
import hotel
h = hotel.hotel([n_rooms]*n_days, P_nom, P_elastic)
def objective(price, hotel):
return -hotel.profit(price)
def constraint(price, hotel): # <= 0
return -min(hotel.remaining(price))
bounds = [P_bounds]*n_days
guess = [P_nom]*n_days
import mystic as my
@my.penalty.quadratic_inequality(constraint, kwds=dict(hotel=h))
def penalty(x):
return 0.0
# try again using a global optimizer
solver = my.solvers.diffev
mon = my.monitors.VerboseMonitor(100)
kwds = dict(disp=True, full_output=True, itermon=mon, npop=40,
args=(h,), gtol=250, ftol=1e-8, maxfun=30000, maxiter=2000)
result = solver(objective, bounds, bounds=bounds, penalty=penalty, **kwds)
print([round(i,2) for i in result[0]])
结果：
>$ python global.py
Generation 0 has Chi-Squared: 3684702.124765
Generation 100 has Chi-Squared: -36493.709890
Generation 200 has Chi-Squared: -36650.498677
Generation 300 has Chi-Squared: -36651.722841
Generation 400 has Chi-Squared: -36651.733274
Generation 500 has Chi-Squared: -36651.733322
Generation 600 has Chi-Squared: -36651.733361
Generation 700 has Chi-Squared: -36651.733361
Generation 800 has Chi-Squared: -36651.733361
STOP("ChangeOverGeneration with {'tolerance': 1e-08, 'generations': 250}")
Optimization terminated successfully.
Current function value: -36651.733361
Iterations: 896
Function evaluations: 24237
[861.07, 893.88, 398.68, 471.4, 9.44, 0.0, 244.67]
我认为为了产生更合理的定价，我会将 P_bounds 值更改为更合理的值。
我正在尝试在 Dynamic room pricing model for hotel revenue management systems 上构建白皮书的实现。以防这个 link 将来死亡，我将粘贴到此处的相关部分：
到目前为止，我目前的实现非常糟糕，因为我真的没有完全理解如何求解非线性最大化方程。
# magical lookup table that returns demand based on those inputs
# this will eventually be a db lookup against past years rental activity and not hardcoded to a specific value
def demand(dateFirstNight, duration):
return 1
# magical function that fetches the price we have allocated for a room on this date to existing customers
# this should be a db lookup against previous stays, and not hardcoded to a specific value
def getPrice(date):
return 75
# Typical room base price
# Defined as: Nominal price of the hotel (usually the average historical price)
nominalPrice = 89
# from the white paper, but perhaps needs to be adjusted in the future using the methods they explain
priceElasticity = 2
# this is an adjustable constant it depends how far forward we want to look into the future when optimizing the prices
# likely this will effect how long this will take to run, so it will be a balancing game with regards to accuracy vs runtime
numberOfDays = 30
def roomsAlocated(dateFirstNight, duration)
roomPriceSum = 0.0
for date in range(dateFirstNight, dateFirstNight+duration-1):
roomPriceSum += getPrice(date)
return demand(dateFirstNight, duration) * (roomPriceSum/(nominalPrice*duration))**priceElasticity
def roomsReserved(date):
# find all stays that contain this date, this
def maximizeRevenue(dateFirstNight):
# we are inverting the price sum which is to be maximized because mystic only does minimization
# and when you minimize the inverse you are maximizing!
return (sum([getPrice(date)*roomsReserved(date) for date in range(dateFirstNight, dateFirstNight+numberOfDays)]))**-1
def constraint(x): # Ol - totalNumberOfRoomsInHotel <= 0
return roomsReserved(date) - totalNumberOfRoomsInHotel
from mystic.penalty import quadratic_inequality
@quadratic_inequality(constraint, k=1e4)
def penalty(x):
return 0.0
from mystic.solvers import diffev2
from mystic.monitors import VerboseMonitor
mon = VerboseMonitor(10)
bounds = [0,1e4]*numberOfDays
result = diffev2(maximizeRevenue, x0=bounds, penalty=penalty, npop=10, gtol=200, disp=False, full_output=True, itermon=mon, maxiter=M*N*100)
任何熟悉使用 mystic 的人都可以给我一些关于如何实现它的建议吗？
虽然您要求使用库 mystic，但在开始非线性优化时，您可能不需要这种细粒度的控制。模块 scipy 应该足够了。以下是一个或多或少的完整解决方案，更正了我认为原始白皮书中关于定价范围的错字：
import numpy as np
from scipy.optimize import minimize
P_nom = 89
P_max = None
price_elasticity = 2
number_of_days = 7
demand = lambda a, L: 1./L
total_rooms = [5]*number_of_days
def objective(P, *args):
return -np.dot(P, O(P, *args))
def worst_leftover(P, C, *args):
return min(np.subtract(C, O(P, *args)))
def X(P, a, L, d, e, P_nom):
return d(a, L)*(sum(P[a:a+L])/P_nom/L)**e
def d(a, L):
return 1.
def O_l(P, l, l_max, *args):
return sum([X(P, a, L, *args)
for a in xrange(0, l)
for L in xrange(l-a+1, l_max+1)])
def O(P, *args):
return [O_l(P, l, *args) for l in xrange(len(P))]
result = minimize(
objective,
[P_nom]*number_of_days,
args=(number_of_days-1, demand, price_elasticity, P_nom),
method='SLSQP',
bounds=[(0, P_max)]*number_of_days,
constraints={
'type': 'ineq',
'fun': worst_leftover,
'args': (total_rooms, number_of_days-1, demand, price_elasticity, P_nom)
},
tol=1e-1,
options={'maxiter': 10**3}
)
print result.x
有几点值得一提：
objective 函数添加了一个减号，用于 scipy 的
minimize()例程，与原始白皮书中引用的最大化形成对比。这将导致result.fun为负数而不是表示总收入。公式好像对参数有点敏感。最小化是正确的（至少，当它说它正确执行时它是正确的——检查
result.success），但如果输入与现实相差太远，那么您可能会发现价格比预期高得多。您可能还想使用比我在以下示例中使用的天数更多的天数。您的白皮书似乎会产生周期性影响。我不太喜欢白皮书的命名方案，因为它与可读代码有关。我改变了一些东西，但有些东西确实很糟糕，应该被替换，比如小写的 l 很容易与数字 1 混淆。
我确实设置了界限，使价格为正而不是负。凭借您的领域专业知识，您应该验证这是正确的决定。
您可能更喜欢比我指定的更严格的公差。这在某种程度上取决于您希望运行时是什么。随意使用
tol参数。此外，对于更严格的公差，您可能会发现options参数中的'maxiter'必须增加才能使minimize()正确收敛。我很确定
total_rooms应该是酒店中尚未预订的房间数量，因为白皮书将其索引为字母 l 而不是常量你有你的原始代码。我将其设置为用于测试目的的常量列表。方法必须是 'SLSQP' 来处理价格范围和房间数量范围。注意不要改这个。
O_l()的计算效率非常低下。如果运行时是个问题，我要采取的第一步是弄清楚如何 cache/memoize 调用X()。所有这些实际上只是第一步，概念验证。它应该没有错误且正确，但它几乎是直接从白皮书中提取的，可以进行一些重构。
任何人，玩得开心，如有任何其他问题，请随时 comment/PM/etc。
抱歉，我来晚了，但我认为接受的答案并没有解决完整的问题，而且还错误地解决了问题。请注意，在局部最小化中，求解接近名义价格的问题并不能给出最佳解决方案。
我们先建一个hotel class:
"""
This file is 'hotel.py'
"""
import math
class hotel(object):
def __init__(self, rooms, price_ave, price_elastic):
self.rooms = rooms
self.price_ave = price_ave
self.price_elastic = price_elastic
def profit(self, P):
# assert len(P) == len(self.rooms)
return sum(j * self._reserved(P, i) for i,j in enumerate(P))
def remaining(self, P): # >= 0
C = self.rooms
# assert len(P) == C
return [C[i] - self._reserved(P, i) for i,j in enumerate(P)]
def _reserved(self, P, day):
max_days = len(self.rooms)
As = range(0, day)
return sum(self._allocated(P, a, L) for a in As
for L in range(day-a+1, max_days+1))
def _allocated(self, P, a, L):
P_nom = self.price_ave
e = self.price_elastic
return math.ceil(self._demand(a, L)*(sum(P[a:a+L])/(P_nom*L))**e)
def _demand(self, a,L): #XXX: fictional non-constant demand function
return abs(1-a)/L + 2*(a**2)/L**2
这是使用 mystic 解决问题的一种方法：
"""
This file is 'local.py'
"""
n_days = 7
n_rooms = 50
P_nom = 85
P_bounds = 0,None
P_elastic = 2
import hotel
h = hotel.hotel([n_rooms]*n_days, P_nom, P_elastic)
def objective(price, hotel):
return -hotel.profit(price)
def constraint(price, hotel): # <= 0
return -min(hotel.remaining(price))
bounds = [P_bounds]*n_days
guess = [P_nom]*n_days
import mystic as my
@my.penalty.quadratic_inequality(constraint, kwds=dict(hotel=h))
def penalty(x):
return 0.0
# using a local optimizer, starting from the nominal price
solver = my.solvers.fmin
mon = my.monitors.VerboseMonitor(100)
kwds = dict(disp=True, full_output=True, itermon=mon,
args=(h,), xtol=1e-8, ftol=1e-8, maxfun=10000, maxiter=2000)
result = solver(objective, guess, bounds=bounds, penalty=penalty, **kwds)
print([round(i,2) for i in result[0]])
结果：
>$ python local.py
Generation 0 has Chi-Squared: -4930.000000
Generation 100 has Chi-Squared: -22353.444547
Generation 200 has Chi-Squared: -22410.402420
Generation 300 has Chi-Squared: -22411.780268
Generation 400 has Chi-Squared: -22413.908944
Generation 500 has Chi-Squared: -22477.869093
Generation 600 has Chi-Squared: -22480.144244
Generation 700 has Chi-Squared: -22480.280379
Generation 800 has Chi-Squared: -22485.563188
Generation 900 has Chi-Squared: -22485.564265
Generation 1000 has Chi-Squared: -22489.341602
Generation 1100 has Chi-Squared: -22489.345912
Generation 1200 has Chi-Squared: -22489.351219
Generation 1300 has Chi-Squared: -22491.994305
Generation 1400 has Chi-Squared: -22491.994518
Generation 1500 has Chi-Squared: -22492.061127
Generation 1600 has Chi-Squared: -22492.573672
Generation 1700 has Chi-Squared: -22492.573690
Generation 1800 has Chi-Squared: -22492.622064
Generation 1900 has Chi-Squared: -22492.622230
Optimization terminated successfully.
Current function value: -22492.622230
Iterations: 1926
Function evaluations: 3346
STOP("CandidateRelativeTolerance with {'xtol': 1e-08, 'ftol': 1e-08}")
[1.15, 20.42, 20.7, 248.1, 220.56, 41.4, 160.09]
再次使用全局优化器：
"""
This file is 'global.py'
"""
n_days = 7
n_rooms = 50
P_nom = 85
P_bounds = 0,None
P_elastic = 2
import hotel
h = hotel.hotel([n_rooms]*n_days, P_nom, P_elastic)
def objective(price, hotel):
return -hotel.profit(price)
def constraint(price, hotel): # <= 0
return -min(hotel.remaining(price))
bounds = [P_bounds]*n_days
guess = [P_nom]*n_days
import mystic as my
@my.penalty.quadratic_inequality(constraint, kwds=dict(hotel=h))
def penalty(x):
return 0.0
# try again using a global optimizer
solver = my.solvers.diffev
mon = my.monitors.VerboseMonitor(100)
kwds = dict(disp=True, full_output=True, itermon=mon, npop=40,
args=(h,), gtol=250, ftol=1e-8, maxfun=30000, maxiter=2000)
result = solver(objective, bounds, bounds=bounds, penalty=penalty, **kwds)
print([round(i,2) for i in result[0]])
结果：
>$ python global.py
Generation 0 has Chi-Squared: 3684702.124765
Generation 100 has Chi-Squared: -36493.709890
Generation 200 has Chi-Squared: -36650.498677
Generation 300 has Chi-Squared: -36651.722841
Generation 400 has Chi-Squared: -36651.733274
Generation 500 has Chi-Squared: -36651.733322
Generation 600 has Chi-Squared: -36651.733361
Generation 700 has Chi-Squared: -36651.733361
Generation 800 has Chi-Squared: -36651.733361
STOP("ChangeOverGeneration with {'tolerance': 1e-08, 'generations': 250}")
Optimization terminated successfully.
Current function value: -36651.733361
Iterations: 896
Function evaluations: 24237
[861.07, 893.88, 398.68, 471.4, 9.44, 0.0, 244.67]
我认为为了产生更合理的定价，我会将 P_bounds 值更改为更合理的值。