00001
00016 #include <iomanip>
00017 #include "flopc.hpp"
00018
00019
00020 #include "OsiClpSolverInterface.hpp"
00021 #define OSI_SOLVER_INTERFACE OsiClpSolverInterface
00022
00023 using namespace std;
00024 using namespace flopc;
00025
00026
00028
00033 class StageNode {
00034 public:
00035 MP_set ASSETS;
00036 MP_index a;
00037 MP_variable x;
00038 MP_expression objFunction;
00039
00040 StageNode *ptParent;
00041 vector<StageNode *> children;
00042 int nodeNmb;
00043 double prob;
00044
00045 StageNode(StageNode *ptPred, const int nmbAssets, const int nodeN,
00046 const double condProb)
00047 : ASSETS(nmbAssets), x(ASSETS), ptParent(ptPred), nodeNmb(nodeN),
00048 prob(condProb)
00049 {
00050 if (ptParent != NULL) {
00051 prob *= ptParent->prob;
00052 ptParent->children.push_back(this);
00053 }
00054 }
00055 virtual ~StageNode(){}
00056
00057 protected:
00059
00061 virtual void make_obj_function_() {
00062 assert (children.size() > 0 && "Leaves should never call this function");
00063
00064 for (int i=0; i < (int) children.size(); i++) {
00065 children[i]->make_obj_function_();
00066 }
00067
00068
00069 objFunction = children[0]->objFunction;
00070 for (int i=1; i < (int) children.size(); i++) {
00071 objFunction = objFunction + children[i]->objFunction;
00072 }
00073 }
00074 };
00075
00077 class RootNode : public StageNode {
00078 public:
00079 MP_constraint initialBudget;
00080
00081 RootNode(const int nmbAssets, const double initWealth)
00082 : StageNode(NULL, nmbAssets, 0, 1.0)
00083 {
00084 initialBudget() = sum(ASSETS, x(ASSETS)) == initWealth;
00085 }
00086
00088
00090 void make_objective_function() {
00091 make_obj_function_();
00092 }
00093 };
00094
00096 class MidStageNode : public StageNode {
00097 public:
00098 MP_constraint cashFlowBalance;
00099 MP_data Return;
00100
00102
00108 MidStageNode(StageNode *ptPred, const int nodeN, const double condProb,
00109 double *ptRetVect)
00110 : StageNode(ptPred, ptPred->ASSETS.size(), nodeN, condProb),
00111 Return(ptRetVect, ASSETS)
00112 {
00113
00114 cashFlowBalance = sum(ASSETS(a), ptParent->x(a) * Return(a))
00115 == sum(ASSETS(a), x(a));
00116 }
00117 virtual ~MidStageNode(){}
00118 };
00119
00121 class LeafNode : public StageNode {
00122 public:
00123 MP_variable w;
00124 MP_variable y;
00125 MP_constraint finalBalance;
00126 MP_data Return;
00127
00129
00133 LeafNode(StageNode *ptPred, const int nodeN, const double condProb,
00134 const double *ptRetVect, const double capTarget)
00135 : StageNode(ptPred, ptPred->ASSETS.size(), nodeN, condProb), Return(ASSETS)
00136 {
00137 Return.value(ptRetVect);
00138
00139 finalBalance = sum(ASSETS, ptParent->x(ASSETS) * Return(ASSETS))
00140 + w() - y() == capTarget;
00141 }
00142
00143 protected:
00145 void make_obj_function_() {
00146 objFunction = prob * (1.3 * w() - 1.1 * y());
00147 }
00148 };
00149
00150
00151
00152
00154 class ScenTreeStruct {
00155 public:
00156 int nmbNodes;
00157 int firstLeaf;
00158
00160 ScenTreeStruct(const int nNodes, const int firstL)
00161 : nmbNodes(nNodes), firstLeaf(firstL)
00162 {}
00163 virtual ~ScenTreeStruct(){}
00164
00166
00169 virtual int get_parent(int n) const = 0;
00170 };
00171
00173 class BinTreeStruct : public ScenTreeStruct {
00174 public:
00176 BinTreeStruct(const int nmbStages)
00177 : ScenTreeStruct(pow(2.0, nmbStages) - 1, pow(2.0, nmbStages-1) - 1) {}
00178
00179 int get_parent(int n) const {
00180 return (n-1) / 2;
00181 }
00182 };
00183
00184
00185
00186 int main()
00187 {
00188
00189
00190 BinTreeStruct scTree(4);
00191
00192
00193 enum {stocks, bonds, nmbAssets};
00194 double InitBudget = 55;
00195 double CapTarget = 80;
00196
00197
00198 double retData[][nmbAssets] = {
00199 {1.25, 1.14},
00200 {1.06, 1.16},
00201 {1.21, 1.17},
00202 {1.07, 1.12},
00203 {1.15, 1.18},
00204 {1.06, 1.12},
00205 {1.26, 1.13},
00206 {1.07, 1.14},
00207 {1.25, 1.15},
00208 {1.06, 1.12},
00209 {1.05, 1.17},
00210 {1.06, 1.15},
00211 {1.05, 1.14},
00212 {1.06, 1.12}
00213 };
00214
00215
00216
00217 MP_model &model = MP_model::getDefaultModel();
00218 model.setSolver(new OSI_SOLVER_INTERFACE);
00219 #ifdef NDEBUG
00220 model.silent();
00221 #else
00222 model.verbose();
00223 #endif
00224
00225
00226 vector<StageNode *> treeNodes;
00227 treeNodes.push_back(new RootNode(nmbAssets, InitBudget));
00228 for (int i = 1; i < scTree.nmbNodes; i++) {
00229 int pred = scTree.get_parent(i);
00230 if (i < scTree.firstLeaf) {
00231
00232 treeNodes.push_back(
00233
00234
00235 new MidStageNode(treeNodes[pred], i, 1.0 / 2, retData[i-1])
00236 );
00237 } else {
00238
00239 treeNodes.push_back(
00240
00241 new LeafNode(treeNodes[pred], i, 1.0 / 2, retData[i-1], CapTarget)
00242 );
00243 }
00244 }
00245
00246 #ifndef NDEBUG
00247
00248 for (int n = 0; n < scTree.nmbNodes; n++) {
00249 cout << "Node " << setw(2) << treeNodes[n]->nodeNmb << " has "
00250 << "prob = " << setw(5) << treeNodes[n]->prob << " and "
00251 << treeNodes[n]->children.size() << " children: ";
00252 for (int cn = 0; cn < (int) treeNodes[n]->children.size(); cn++)
00253 cout << setw(2) << treeNodes[n]->children[cn]->nodeNmb << " ";
00254 cout << endl;
00255 }
00256 #endif
00257
00258
00259 assert (dynamic_cast<RootNode *>(treeNodes[0]) != NULL
00260 && "The first node in treeNodes should be the Root");
00261 RootNode *ptRoot = dynamic_cast<RootNode *>(treeNodes[0]);
00262
00263 ptRoot->make_objective_function();
00264 model.setObjective(ptRoot->objFunction);
00265 model.attach();
00266 model->writeMps("investment");
00267
00268
00269 model.minimize();
00270
00271
00272 cout << endl << "Optimal objective value = " << model->getObjValue() << endl
00273 << "First-stage x:";
00274 ptRoot->x.display();
00275
00276 for (int n = 0; n < (int) treeNodes.size(); n++)
00277 delete treeNodes[n];
00278
00279 return 0;
00280 }