2245.sagews

"""
Statement of Problem:

Miranda and Fackler

Dairy cow properties:
 N1 lacation cycles:      #s in set S {1,2,N1}
 N2 productivity classes  #x in set X {1,2,N2}
 Each cow belonging to productivity class x yields Qx * y(s) tons of milk during the Sth lactation cycle
 The farmer does not know the productivity class of a cow until after its first lacation

 State Variables:
  Dairy cows lactation cycle:  If Zt = 0 then St+1 = St + 1 else St+1 = 1;   
  Stochastic quality of cow:   If Zt = 0 then Xt+1 = Xt with probability 1 else Xt+1 = Xi with probability Wi where Wi is prob. 
                               of getting a cow ofclass i.   
 
 Choice Variables:
  z = 0 (keep cow); z = 1 (replace cow);

 Benefit Function: 
  Profit(z,x,s) = { If z = 0 p * Qx * y(s) else p * Qx * y(s) - c} c is cost of replacing a cow... note: replace cow after you milk it...

 Formal Statement of Problem: 
  max(Zt = 0,1) sum from t= 0 to infinity [Profit(Zt,XtSt)/(1+r)^t] such that: 
    If Zt = 0 then Xt+1 = Xt with probability 1 else Xt+1 = Xi with probability Wi where Wi is prob. of getting a cow ofclass i. 
    If Zt = 0 then St+1 = St + 1 else St+1 = 1; 

 Bellman's Equation:
  V(St,Xt) = max(Zt = 0,1) {If Zt = 0 then p * Qx * y(st) + B*V(s+1,Xt) else p * Qx * y(st) - c + B * (sum from i=1 to N2 (Wi * V(1,Xi))}
  We will let the Qx * y(St) just be a logarithmic function for now. 
"""

import numpy; #gives access to stuff like arrays, zeros, etc... 
import scipy;

#class to calculate a value function iteration 
class ValueFunctionIteration(): #define class

 #begin constructor
 def __init__(self,maxIterations):

  self.converged= false;
  self.currentIteration = 0;
  self.delta= .9;
  self.maxIterations = maxIterations;
  self.N1= 10;
  self.N2= 10;
  self.c = 10;
  self.p = 5; # price of milk is $5
  self.tau = .3;

  #set VRHS(x,s) = 0 for every state x in X, s in S;
  self.valueRightHandSide = numpy.matrix(numpy.zeros([self.N1,self.N2],float));
  self.valueLeftHandSide = numpy.matrix(numpy.zeros([self.N1,self.N2],float));
  
  self.policyFunction =[];
  self.vKeepCow = 0;
  self.vReplaceCow= 0;

  #probabilities of getting a cow of a certain quality
  self.cowQualityProbabilities=[];
  for i in range(0,self.N2):        
      self.cowQualityProbabilities.append(1/self.N2);

  print 'cow quality proabilities:', self.cowQualityProbabilities;

 #end constructor

 #begin functions

 def performCalculation(self):
 
    while (self.currentIteration < self.maxIterations):  
             
         #for every quality state Xt in X 
         for x in range(0,self.N1 - 1):
           # for every lactation age St in S
           for s in range(0,self.N2 - 1):
         
             #value function value if you keep the cow...
             vKeepCow = float(self.p * log(s) + self.delta * self.valueRightHandSide[s + 1,x]);

             #calculate expected value of getting a new cow...
             expectedValue= 0;
             for k in range(0,len(self.cowQualityProbabilities)):
               expectedValue += k * self.valueRightHandSide[1,x];
               #print expectedValue;

             #value function if you replace the cow...
             vReplaceCow= float(self.p * log(s) - self.c + self.delta * expectedValue);

             if vReplaceCow > vKeepCow:
               self.policyFunction.append('replace cow');
               self.valueLeftHandSide[x,s] = vReplaceCow;              
             else:
               self.policyFunction.append('keep cow');
               self.valueLeftHandSide[x,s] = vKeepCow;             
                                      
             s += 1;             
           x += 1; 

         difference= numpy.linalg.norm(self.valueLeftHandSide - self.valueRightHandSide);         
         print difference;

         if difference < self.tau:
           self.converged= true;
           print 'converged';
           break;
         else: 
           self.valueRightHandSide= self.valueLeftHandSide;
           self.currentIteration += 1;         
    
    print 'converged: ', self.converged;  
     
 #end functions

#create a class instance and execute problem
test= ValueFunctionIteration(20);
test.performCalculation();

cow quality proabilities: [1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10] inf nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan converged: False