Maximize sum subject to conditions in Python mystic

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?

Solution

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.