RMOL Logo Get Revenue Management Optimisation Library at SourceForge.net. Fast, secure and Free Open Source software downloads

DPOptimiser.cpp

Go to the documentation of this file.
00001 // //////////////////////////////////////////////////////////////////////
00002 // Import section
00003 // //////////////////////////////////////////////////////////////////////
00004 // GSL Random Number Generation (GSL Reference Manual, version 1.7, Chapter 19)
00005 #include <gsl/gsl_cdf.h>
00006 #include <gsl/gsl_randist.h>
00007 // C
00008 #include <assert.h>
00009 // STL
00010 #include <sstream>
00011 #include <vector>
00012 #include <cmath>
00013 // RMOL
00014 #include <rmol/basic/BasConst_General.hpp>
00015 #include <rmol/bom/DPOptimiser.hpp>
00016 #include <rmol/bom/Bucket.hpp>
00017 #include <rmol/bom/BucketHolder.hpp>
00018 #include <rmol/service/Logger.hpp>
00019 
00020 namespace RMOL {
00021   
00022   // ////////////////////////////////////////////////////////////////////
00023   void DPOptimiser::
00024   optimalOptimisationByDP (const ResourceCapacity_T iCabinCapacity,
00025                            BucketHolder& ioBucketHolder,
00026                            BidPriceVector_T& ioBidPriceVector) {
00027      // Number of classes/buckets: n
00028     const short nbOfClasses = ioBucketHolder.getSize();
00029 
00030     // Number of values of x to compute for each Vj(x).
00031     const int maxValue = static_cast<int> (iCabinCapacity * DEFAULT_PRECISION);
00032 
00033     // Vector of the Expected Maximal Revenue (Vj).
00034     std::vector< std::vector<double> > MERVectorHolder;
00035 
00036     // Vector of V_0(x).
00037     std::vector<double> initialMERVector (maxValue+1, 0.0);
00038     MERVectorHolder.push_back (initialMERVector);
00039 
00040     // Current cumulative protection level (y_j * DEFAULT_PRECISION).
00041     // Initialise with y_0 = 0.
00042     int currentProtection = 0;
00043 
00044     int currentBucketIndex = 1;
00045     ioBucketHolder.begin();
00046     
00047     while (currentProtection < maxValue && currentBucketIndex < nbOfClasses) {
00048     //while (currentBucketIndex == 1) {
00049       bool protectionChanged = false;
00050       double nextProtection = 0.0;
00051       std::vector<double> currentMERVector;
00052       // double testGradient = 10000;
00053       
00054       Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00055       const double meanDemand = currentBucket.getMean();
00056       const double SDDemand = currentBucket.getStandardDeviation();
00057       const double currentYield = currentBucket.getAverageYield();
00058       const double errorFactor = 1;//gsl_cdf_gaussian_Q (-meanDemand, SDDemand);
00059 
00060       Bucket& nextBucket = ioBucketHolder.getNextBucket();
00061       const double nextYield = nextBucket.getAverageYield();
00062       
00063       // For x <= currentProtection (y_(j-1)), V_j(x) = V_(j-1)(x).
00064       for (int x = 0; x <= currentProtection; ++x) {
00065         const double MERValue = MERVectorHolder.at(currentBucketIndex-1).at(x);
00066         currentMERVector.push_back(MERValue);
00067       }
00068       
00069       // Vector of gaussian pdf values.
00070       std::vector<double> pdfVector;
00071       for (int s = 0; s <= maxValue - currentProtection; ++s) {
00072         const double pdfValue =
00073           gsl_ran_gaussian_pdf (s/DEFAULT_PRECISION - meanDemand, SDDemand);
00074         pdfVector.push_back(pdfValue);
00075       }
00076 
00077       // Vector of gaussian cdf values.
00078       std::vector<double> cdfVector;
00079       for (int s = 0; s <= maxValue - currentProtection; ++s) {
00080         const double cdfValue =
00081           cdfGaussianQ (s/DEFAULT_PRECISION - meanDemand, SDDemand);
00082         cdfVector.push_back(cdfValue);
00083       }
00084       
00085       // Compute V_j(x) for x > currentProtection (y_(j-1)).
00086       for (int x = currentProtection + 1; x <= maxValue; ++x) {
00087         const double lowerBound = static_cast<double> (x - currentProtection);
00088         
00089         // Compute the first integral in the V_j(x) formulation (see
00090         // the memo of Jerome Contant).
00091         const double power1 = - 0.5 * meanDemand * meanDemand /
00092           (SDDemand * SDDemand);
00093         const double e1 = std::exp (power1);
00094         const double power2 = 
00095           - 0.5 * (lowerBound / DEFAULT_PRECISION - meanDemand) *
00096           (lowerBound / DEFAULT_PRECISION - meanDemand) /
00097           (SDDemand * SDDemand);
00098         const double e2 = exp (power2);
00099         /*
00100         const double integralResult1 = currentYield * 
00101           ((e1 - e2) * SDDemand / sqrt (2 * 3.14159265) +
00102            meanDemand * gsl_cdf_gaussian_Q (-meanDemand, SDDemand) -
00103            meanDemand * gsl_cdf_gaussian_Q (lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand));
00104         */
00105         const double integralResult1 = currentYield * 
00106           ((e1 - e2) * SDDemand / sqrt (2 * 3.14159265) +
00107            meanDemand * cdfGaussianQ (-meanDemand, SDDemand) -
00108            meanDemand * cdfGaussianQ (lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand));
00109         
00110         double integralResult2 = 0.0;
00111         
00112         for (int s = 0; s < lowerBound; ++s) {
00113           const double partialResult =
00114             2 * MERVectorHolder.at(currentBucketIndex-1).at(x-s) *
00115             pdfVector.at(s);
00116           
00117           integralResult2 += partialResult;
00118         }
00119         integralResult2 -= MERVectorHolder.at(currentBucketIndex-1).at(x) *
00120           pdfVector.at(0);
00121         
00122         const int intLowerBound = static_cast<int>(lowerBound);
00123         integralResult2 += 
00124           MERVectorHolder.at(currentBucketIndex-1).at(x - intLowerBound) *
00125           pdfVector.at(intLowerBound);
00126         
00127         integralResult2 /= 2 * DEFAULT_PRECISION;
00128         /*
00129         for (int s = 0; s < lowerBound; ++s) {
00130           const double partialResult =
00131             (MERVectorHolder.at(currentBucketIndex-1).at(x-s) +
00132              MERVectorHolder.at(currentBucketIndex-1).at(x-s-1)) *
00133             (cdfVector.at(s+1) - cdfVector.at(s)) / 2;
00134           integralResult2 += partialResult;
00135         }
00136         */
00137         const double firstElement = integralResult1 + integralResult2;
00138         
00139         // Compute the second integral in the V_j(x) formulation (see
00140         // the memo of Jerome Contant).
00141         const double constCoefOfSecondElement =
00142           currentYield * lowerBound / DEFAULT_PRECISION +
00143           MERVectorHolder.at(currentBucketIndex-1).at(currentProtection);
00144         const double secondElement = constCoefOfSecondElement * 
00145     //gsl_cdf_gaussian_Q(lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand);
00146           cdfGaussianQ (lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand);
00147         const double MERValue = (firstElement + secondElement) / errorFactor;
00148 
00149         
00150         assert (currentMERVector.size() > 0);
00151         const double lastMERValue = currentMERVector.back();
00152 
00153         const double currentGradient = 
00154           (MERValue - lastMERValue) * DEFAULT_PRECISION;
00155 
00156         //assert (currentGradient >= 0);
00157         if (currentGradient < -0) {
00158           std::ostringstream ostr;
00159           ostr << currentGradient  << std::endl
00160                << "x = " << x << std::endl
00161                << "class: " << currentBucketIndex << std::endl;
00162           RMOL_LOG_DEBUG (ostr.str());
00163         }
00164           
00165         /*
00166         assert (currentGradient <= testGradient);
00167         testGradient = currentGradient;
00168         */
00169         if (protectionChanged == false && currentGradient <= nextYield) {
00170           nextProtection = x - 1;
00171           protectionChanged = true;
00172         }
00173 
00174          if (protectionChanged == true && currentGradient > nextYield) {
00175           protectionChanged = false;
00176         }
00177         
00178         if (protectionChanged == false && x == maxValue) {
00179           nextProtection = maxValue;
00180         }
00181         
00182         currentMERVector.push_back (MERValue); 
00183       }
00184 
00185       // DEBUG
00186       RMOL_LOG_DEBUG ("Vmaxindex = " << currentMERVector.back());
00187         
00188       MERVectorHolder.push_back (currentMERVector);
00189      
00190       const double realProtection = nextProtection / DEFAULT_PRECISION;
00191       const double bookingLimit = iCabinCapacity - realProtection;
00192 
00193       currentBucket.setCumulatedProtection (realProtection);
00194       nextBucket.setCumulatedBookingLimit (bookingLimit);
00195 
00196       currentProtection = static_cast<int> (std::floor (nextProtection));
00197       
00198       ioBucketHolder.iterate();
00199       ++currentBucketIndex;
00200     }
00201   }
00202 
00203   // ////////////////////////////////////////////////////////////////////
00204   double DPOptimiser::cdfGaussianQ (const double c, const double sd) {
00205     const double power = - c * c * 0.625 / (sd * sd);
00206     const double e = std::sqrt (1 - std::exp (power));
00207     double result = 0.0;
00208     
00209     if (c >= 0) {
00210       result = 0.5 * (1 - e);
00211 
00212     } else {
00213       result = 0.5 * (1 + e);
00214     }
00215     
00216     return result;
00217   }
00218   
00219 }
SourceForge Logo

Generated on Wed Feb 9 2011 17:09:14 for RMOL by Doxygen 1.7.1