2144.sagews

"""
Class to calculate a bellmans discrete time dynamic optimization problem
"""
class BellmansDiscreteTime():    # define class
 
 # Begin Constructor 

 def __init__(self,n):     

  self.n = n  # number of iterations

  # I used these because I couldn't do calculus on indexed variables, so its a little bit of a workaround
  self.xt = var('xt') 
  self.yt = var('yt')

  # Starting objective function... Will modify this as the periods go on...
  # Note: ** is how you do exponents in python
  # Function looks like: 10xt - .1yt^2 where t is a subscript denoting period
  self.origObjectiveF = (10 * self.xt) + (-.1 * self.yt**2) #initial objective function
  self.objectiveF = self.origObjectiveF
   
 # End Constructor

 # Begin Functions

 # Python doesn't support overloaded functions, so... you have original and origional + 1
 def performCalculation(self):
  return self.performCalculation1(self.n - 1)

 def performCalculation1(self,r):

  print self.objectiveF
 
  f(xt,yt) = self.objectiveF
  Dyt = f.derivative(yt)
  ytSolved= solve(Dyt,yt,solution_dict=True)        

  myYt = var('myYt')
  for soln in ytSolved: myYt = soln[yt]

  if r == 0:
   
   # xt = 0 because resource is used up
   # basically here we are on last period
   # we are just doing the final substitution and closing up the recursion

   f= f.substitute(xt = 0,yt = myYt)
   print f
   print 'done with recursion'
       
  else:        

   # here we are performing the substitutions and building the objective function
   # that we will be starting with in the next period... 
   # finally we call ourself to perform the recursion
 
   f = f.substitute(yt = myYt)   
   f= f.substitute(xt= (xt + yt))
   self.objectiveF = self.origObjectiveF + f    
   self.performCalculation1(r-1)        
  
  return f
 
 # End Functions

# Create a class instance, and execute problem
test= BellmansDiscreteTime(5) # create class instance with 5 periods... n = 5;
print "n= %d" % (test.n, )  
test.performCalculation()

n= 5 -0.100000000000000*yt^2 + 10*xt (xt, yt) |--> -0.100000000000000*yt^2 + 20*xt + 10*yt (xt, yt) |--> -0.100000000000000*yt^2 + 30*xt + 20*yt + 250.000000000000 (xt, yt) |--> -0.100000000000000*yt^2 + 40*xt + 30*yt + 1250.00000000000 (xt, yt) |--> -0.100000000000000*yt^2 + 50*xt + 40*yt + 3500.00000000000 (xt, yt) |--> 7500.00000000000 done with recursion

\newcommand{\Bold}[1]{\mathbf{#1}}\left( \mbox{xt}, \mbox{yt} \right) \ {\mapsto} \ 10 \, \mbox{xt} + 10 \, \mbox{yt}