RMOL Logo  1.00.9
C++ library of Revenue Management and Optimisation classes and functions
Loading...
Searching...
No Matches
DPOptimiser.cpp
Go to the documentation of this file.
1// //////////////////////////////////////////////////////////////////////
2// Import section
3// //////////////////////////////////////////////////////////////////////
4// STL
5#include <cassert>
6#include <sstream>
7#include <vector>
8#include <cmath>
9// Boost Math
10#include <boost/math/distributions/normal.hpp>
11// StdAir
12#include <stdair/bom/LegCabin.hpp>
13#include <stdair/bom/VirtualClassStruct.hpp>
14#include <stdair/service/Logger.hpp>
15// RMOL
18
19namespace RMOL {
20
21 // ////////////////////////////////////////////////////////////////////
22 void DPOptimiser::optimalOptimisationByDP (stdair::LegCabin& ioLegCabin) {
23 // // Number of classes/buckets: n
24 // const short nbOfClasses = ioBucketHolder.getSize();
25
26 // // Number of values of x to compute for each Vj(x).
27 // const int maxValue = static_cast<int> (iCabinCapacity * DEFAULT_PRECISION);
28
29 // // Vector of the Expected Maximal Revenue (Vj).
30 // std::vector< std::vector<double> > MERVectorHolder;
31
32 // // Vector of V_0(x).
33 // std::vector<double> initialMERVector (maxValue+1, 0.0);
34 // MERVectorHolder.push_back (initialMERVector);
35
36 // // Current cumulative protection level (y_j * DEFAULT_PRECISION).
37 // // Initialise with y_0 = 0.
38 // int currentProtection = 0;
39
40 // int currentBucketIndex = 1;
41 // ioBucketHolder.begin();
42
43 // while (currentProtection < maxValue && currentBucketIndex < nbOfClasses) {
44 // //while (currentBucketIndex == 1) {
45 // bool protectionChanged = false;
46 // double nextProtection = 0.0;
47 // std::vector<double> currentMERVector;
48 // // double testGradient = 10000;
49
50 // Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
51 // const double meanDemand = currentBucket.getMean();
52 // const double SDDemand = currentBucket.getStandardDeviation();
53 // const double currentYield = currentBucket.getAverageYield();
54 // const double errorFactor = 1.0;
55
56 // Bucket& nextBucket = ioBucketHolder.getNextBucket();
57 // const double nextYield = nextBucket.getAverageYield();
58
59 // // For x <= currentProtection (y_(j-1)), V_j(x) = V_(j-1)(x).
60 // for (int x = 0; x <= currentProtection; ++x) {
61 // const double MERValue = MERVectorHolder.at(currentBucketIndex-1).at(x);
62 // currentMERVector.push_back (MERValue);
63 // }
64
65 // //
66 // boost::math::normal lNormalDistribution (meanDemand, SDDemand);
67
68 // // Vector of gaussian pdf values.
69 // std::vector<double> pdfVector;
70 // for (int s = 0; s <= maxValue - currentProtection; ++s) {
71 // const double pdfValue =
72 // boost::math::pdf (lNormalDistribution, s/DEFAULT_PRECISION);
73 // pdfVector.push_back (pdfValue);
74 // }
75
76 // // Vector of gaussian cdf values.
77 // std::vector<double> cdfVector;
78 // for (int s = 0; s <= maxValue - currentProtection; ++s) {
79 // const double cdfValue =
80 // boost::math::cdf (boost::math::complement (lNormalDistribution,
81 // s/DEFAULT_PRECISION));
82 // cdfVector.push_back (cdfValue);
83 // }
84
85 // // Compute V_j(x) for x > currentProtection (y_(j-1)).
86 // for (int x = currentProtection + 1; x <= maxValue; ++x) {
87 // const double lowerBound = static_cast<double> (x - currentProtection);
88
89 // // Compute the first integral in the V_j(x) formulation (see
90 // // the memo of Jerome Contant).
91 // const double power1 =
92 // - 0.5 * meanDemand * meanDemand / (SDDemand * SDDemand);
93 // const double e1 = std::exp (power1);
94 // const double power2 =
95 // - 0.5 * (lowerBound / DEFAULT_PRECISION - meanDemand)
96 // * (lowerBound / DEFAULT_PRECISION - meanDemand)
97 // / (SDDemand * SDDemand);
98 // const double e2 = std::exp (power2);
99
100 // const double cdfValue0 =
101 // boost::math::cdf (boost::math::complement (lNormalDistribution,
102 // 0.0));
103 // const double cdfValue1 =
104 // boost::math::cdf(boost::math::complement(lNormalDistribution,
105 // lowerBound/DEFAULT_PRECISION));
106 // const double integralResult1 = currentYield
107 // * ((e1 - e2) * SDDemand / sqrt (2 * 3.14159265)
108 // + meanDemand * (cdfValue0 - cdfValue1));
109
110 // double integralResult2 = 0.0;
111
112 // for (int s = 0; s < lowerBound; ++s) {
113 // const double partialResult =
114 // 2 * MERVectorHolder.at(currentBucketIndex-1).at(x-s)
115 // * pdfVector.at(s);
116
117 // integralResult2 += partialResult;
118 // }
119 // integralResult2 -= MERVectorHolder.at(currentBucketIndex-1).at(x) *
120 // pdfVector.at(0);
121
122 // const int intLowerBound = static_cast<int>(lowerBound);
123 // integralResult2 +=
124 // MERVectorHolder.at(currentBucketIndex-1).at(x - intLowerBound) *
125 // pdfVector.at(intLowerBound);
126
127 // integralResult2 /= 2 * DEFAULT_PRECISION;
128 // /*
129 // for (int s = 0; s < lowerBound; ++s) {
130 // const double partialResult =
131 // (MERVectorHolder.at(currentBucketIndex-1).at(x-s) +
132 // MERVectorHolder.at(currentBucketIndex-1).at(x-s-1)) *
133 // (cdfVector.at(s+1) - cdfVector.at(s)) / 2;
134 // integralResult2 += partialResult;
135 // }
136 // */
137 // const double firstElement = integralResult1 + integralResult2;
138
139 // // Compute the second integral in the V_j(x) formulation (see
140 // // the memo of Jerome Contant).
141 // const double constCoefOfSecondElement =
142 // currentYield * lowerBound / DEFAULT_PRECISION
143 // + MERVectorHolder.at(currentBucketIndex-1).at(currentProtection);
144
145 // const double secondElement = constCoefOfSecondElement
146 // * boost::math::cdf(boost::math::complement(lNormalDistribution,
147 // lowerBound/DEFAULT_PRECISION));
148
149 // const double MERValue = (firstElement + secondElement) / errorFactor;
150
151
152 // assert (currentMERVector.size() > 0);
153 // const double lastMERValue = currentMERVector.back();
154
155 // const double currentGradient =
156 // (MERValue - lastMERValue) * DEFAULT_PRECISION;
157
158 // //assert (currentGradient >= 0);
159 // if (currentGradient < -0) {
160 // std::ostringstream ostr;
161 // ostr << currentGradient << std::endl
162 // << "x = " << x << std::endl
163 // << "class: " << currentBucketIndex << std::endl;
164 // STDAIR_LOG_DEBUG (ostr.str());
165 // }
166
167 // /*
168 // assert (currentGradient <= testGradient);
169 // testGradient = currentGradient;
170 // */
171 // if (protectionChanged == false && currentGradient <= nextYield) {
172 // nextProtection = x - 1;
173 // protectionChanged = true;
174 // }
175
176 // if (protectionChanged == true && currentGradient > nextYield) {
177 // protectionChanged = false;
178 // }
179
180 // if (protectionChanged == false && x == maxValue) {
181 // nextProtection = maxValue;
182 // }
183
184 // currentMERVector.push_back (MERValue);
185 // }
186
187 // // DEBUG
188 // STDAIR_LOG_DEBUG ("Vmaxindex = " << currentMERVector.back());
189
190 // MERVectorHolder.push_back (currentMERVector);
191
192 // const double realProtection = nextProtection / DEFAULT_PRECISION;
193 // const double bookingLimit = iCabinCapacity - realProtection;
194
195 // currentBucket.setCumulatedProtection (realProtection);
196 // nextBucket.setCumulatedBookingLimit (bookingLimit);
197
198 // currentProtection = static_cast<int> (std::floor (nextProtection));
199
200 // ioBucketHolder.iterate();
201 // ++currentBucketIndex;
202 // }
203 }
204
205}
static void optimalOptimisationByDP(stdair::LegCabin &)