Maximize sum subject to conditions in Python mystic

Question

I am trying to build an implementation of a white paper on Dynamic room pricing model for hotel revenue management systems. In case this link dies in the future, i am pasting in the relevant section here:

My current implmentation thus far is quite massively broken, as I really do not fully comprehend how to solve non-linear maximization equations.

# 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)

Can anyone that is familiar with working with mystic give me some advice on how to implement this?

Answer 1

While you asked for use of the library mystic, you probably don't need such fine-grained control when starting out with nonlinear optimization. The module scipy ought to suffice. What follows is a more-or-less complete solution, correcting what I perceive to be a typo in the original whitepaper regarding the pricing bounds:

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

A couple of points worth mentioning:

1. The objective function has an added minus sign for use with scipy's minimize() routine, contrasting with the maximization referenced in the original whitepaper. This will cause result.fun to be negative rather than indicating total revenue.

2. The formula seems to be a bit sensitive to the parameters. The minimization is correct (at least, it's correct when it says it executed correctly -- check result.success), but if the inputs are too far off from reality then you're likely to find prices much higher than anticipated. You also probably want to use more days than I did in the following sample. There seems to be something like a periodic effect your whitepaper induces.

3. I'm not really a fan of the whitepaper's naming scheme as it relates to readable code. I changed a few things, but some things are truly atrocious and ought to be replaced, like a lower-case l which could be easily confused for the number 1.

4. I did set the bounds so that prices are positive rather than negative. With your domain expertise, you should verify that was the correct decision.

5. You might prefer tighter tolerances than I specified. That depends somewhat on what you want runtimes to be. Feel free to play with the tol parameter. Additionally, with tighter tolerances you might find that 'maxiter' in the options parameter has to be increased for minimize() to converge properly.

6. I'm pretty sure total_rooms is supposed to be the number of not-yet-booked rooms in the hotel since the whitepaper has it indexed by the letter l rather than constant as you had in your original code. I set it to be a constant list for testing purposes.

7. The method will have to be 'SLSQP' to deal with the bounds on the prices and the bounds on the number of rooms. Take care not to change this one.

8. There is a massive inefficiency in how O_l() is computed. If runtime is an issue, the first step I would take is figuring out how to cache/memoize calls to X(). All of this is really just a first pass, proof-of-concept. It should be reasonably bug-free and correct, but it is pulled nearly directly from a whitepaper and could do with some re-factoring.

Anywho, have fun and feel free to comment/PM/etc with any further questions.

Answer 2

Sorry I'm late to answer this, but I think the accepted answer is not solving the complete problem, and furthermore solving it incorrectly. Note how in the local minimization, solving near nominal price doesn't give the best solution.

Let's first build a 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

Here's one way to solve it using 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]])

Results:

>$ 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]

Here it is again, using a global optimizer:

"""
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]])

Results:

>$ 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]

I think for this to yield more reasonable pricing, I'd change P_bounds values to something more reasonable.