Go to the documentation of this file.00001
00002
00003
00004
00005 #include <gsl/gsl_cdf.h>
00006 #include <gsl/gsl_randist.h>
00007
00008 #include <assert.h>
00009
00010 #include <sstream>
00011 #include <vector>
00012 #include <cmath>
00013
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
00028 const short nbOfClasses = ioBucketHolder.getSize();
00029
00030
00031 const int maxValue = static_cast<int> (iCabinCapacity * DEFAULT_PRECISION);
00032
00033
00034 std::vector< std::vector<double> > MERVectorHolder;
00035
00036
00037 std::vector<double> initialMERVector (maxValue+1, 0.0);
00038 MERVectorHolder.push_back (initialMERVector);
00039
00040
00041
00042 int currentProtection = 0;
00043
00044 int currentBucketIndex = 1;
00045 ioBucketHolder.begin();
00046
00047 while (currentProtection < maxValue && currentBucketIndex < nbOfClasses) {
00048
00049 bool protectionChanged = false;
00050 double nextProtection = 0.0;
00051 std::vector<double> currentMERVector;
00052
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;
00059
00060 Bucket& nextBucket = ioBucketHolder.getNextBucket();
00061 const double nextYield = nextBucket.getAverageYield();
00062
00063
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
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
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
00086 for (int x = currentProtection + 1; x <= maxValue; ++x) {
00087 const double lowerBound = static_cast<double> (x - currentProtection);
00088
00089
00090
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
00101
00102
00103
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
00130
00131
00132
00133
00134
00135
00136
00137 const double firstElement = integralResult1 + integralResult2;
00138
00139
00140
00141 const double constCoefOfSecondElement =
00142 currentYield * lowerBound / DEFAULT_PRECISION +
00143 MERVectorHolder.at(currentBucketIndex-1).at(currentProtection);
00144 const double secondElement = constCoefOfSecondElement *
00145
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
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
00167
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
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 }