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

Generated on Fri Sep 11 06:31:59 2009 for RMOL by Doxygen 1.5.8