source: trunk/Cbc/examples/link.cpp @ 1565

Last change on this file since 1565 was 121, checked in by forrest, 14 years ago

new example

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.5 KB
Line 
1// Copyright (C) 2005, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7
8#include <cassert>
9#include <iomanip>
10
11
12// For Branch and bound
13#include "OsiSolverInterface.hpp"
14#include "CbcModel.hpp"
15#include "CoinModel.hpp"
16// For Linked Ordered Sets
17#include "CbcBranchLink.hpp"
18#include "OsiClpSolverInterface.hpp"
19
20#include  "CoinTime.hpp"
21
22
23/************************************************************************
24
25This shows how we can define a new branching method to solve problems with
26nonlinearities and discontinuities.
27
28We are going to solve the problem
29
30minimize  10.0*x + y + z - w
31
32where x, y are continuous between 1 and 10, z can take the values 0, 0.1, 0.2 up to 1.0
33and w can take the values 0 or 1.  There is one constraint
34
35w  <= x*z**2 + y*sqrt(z)
36
37One could try to use logarithms to make the problem separable but that is a very
38weak formulation as we want to branch on z directly.  The answer is the concept of linked
39special ordered sets.  The generalization with column generation can be even more powerful
40but here we are limiting z to discrete values to avoid column generation.
41
42The idea is simple:
43
44A linear variable is a convex combination of its lower bound and upper bound!
45If x must lie between 2 and 10 then we can substitute for x  as x == 2.0*xl + 10.0*xu where
46xl + xu == 1.0.  At first this looks cumbersome but if we have xl0, xl1, ... xl10 and corresponding
47xu and yl and yu then we can write:
48
49x == sum 2.0*xl[i] + 10.0* xu[i] where sum xl[i] + xu[i] == 1.0
50and
51x*z**2 == 0.02*xl1 + 0.1*xu1 + 0.08*xl2 + 0.4*xu2 .... + 2.0*xl10 + 10.0*xu10
52
53with similar substitutions for y and y*sqrt(z)
54
55And now the problem is satisfied if w is 0 or 1 and xl[i], xu[i], yl[i] and yu[i] are only
56nonzero for one i.
57
58So this is just like a special ordered set of type 1 but on four sets simultaneously.
59Also note that convexity requirements for any non-linear functions are not needed.
60
61So we need a new branching method to do that - see CbcBranchLink.?pp
62
63We are going to need a CbcBranchLink method to see whether we are satisfied etc and also to
64create another branching object which knows how to fix variables.  We might be able to use an
65existing method for the latter but let us create two methods CbcLink and
66CbcLinkBranchingObject.
67
68For CbcLink we will need the following methods:
69Constructot/Destructor
70infeasibility - returns 0.0 if feasible otherwise some measure of infeasibility
71feasibleRegion - sets bounds to contain current solution
72createBranch - creates a CbcLinkBranchingObject
73
74For CbcLinkBranchingObject we need:
75Constructor/Destructor
76branch - does actual fixing
77print - optional for debug purposes.
78
79The easiest way to do this is to cut and paste from CbcBranchActual to get current
80SOS stuff and then modify that.
81
82************************************************************************/
83
84int main (int argc, const char *argv[])
85{
86
87  OsiClpSolverInterface solver1;
88
89  // Create model
90  CoinModel build;
91  // Keep x,y and z for reporting purposes in rows 0,1,2
92  // Do these
93  double value=1.0;
94  int row=-1;
95  // x
96  row=0;
97  build.addColumn(1,&row,&value,1.0,10.0,10.0);
98  // y
99  row=1;
100  build.addColumn(1,&row,&value,1.0,10.0,1.0);
101  // z
102  row=2;
103  build.addColumn(1,&row,&value,0.0,1.0,1.0);
104  // w
105  row=3;
106  build.addColumn(1,&row,&value,0.0,1.0,-1.0);
107  build.setInteger(3);
108  // Do columns so we know where each is
109  int i;
110  for (i=4;i<4+44;i++)
111    build.setColumnBounds(i,0.0,1.0);
112  // Now do rows
113  // x
114  build.setRowBounds(0,0.0,0.0);
115  for (i=0;i<11;i++) {
116    // xl
117    build.setElement(0,4+4*i,-1.0);
118    // xu
119    build.setElement(0,4+4*i+1,-10.0);
120  }
121  // y
122  build.setRowBounds(1,0.0,0.0);
123  for (i=0;i<11;i++) {
124    // yl
125    build.setElement(1,4+4*i+2,-1.0);
126    // yu
127    build.setElement(1,4+4*i+3,-10.0);
128  }
129  // z - just use x part
130  build.setRowBounds(2,0.0,0.0);
131  for (i=0;i<11;i++) {
132    // xl
133    build.setElement(2,4+4*i,-0.1*i);
134    // xu
135    build.setElement(2,4+4*i+1,-0.1*i);
136  }
137  // w  <= x*z**2 + y* sqrt(z)
138  build.setRowBounds(3,-COIN_DBL_MAX,0.0);
139  for (i=0;i<11;i++) {
140    double value = 0.1*i;
141    // xl * z**2
142    build.setElement(3,4+4*i,-1.0*value*value);
143    // xu * z**2
144    build.setElement(3,4+4*i+1,-10.0*value*value);
145    // yl * sqrt(z)
146    build.setElement(3,4+4*i+2,-1.0*sqrt(value));
147    // yu * sqrt(z)
148    build.setElement(3,4+4*i+3,-10.0*sqrt(value));
149  }
150  // and convexity for x and y
151  // x
152  build.setRowBounds(4,1.0,1.0);
153  for (i=0;i<11;i++) {
154    // xl
155    build.setElement(4,4+4*i,1.0);
156    // xu
157    build.setElement(4,4+4*i+1,1.0);
158  }
159  // y
160  build.setRowBounds(5,1.0,1.0);
161  for (i=0;i<11;i++) {
162    // yl
163    build.setElement(5,4+4*i+2,1.0);
164    // yu
165    build.setElement(5,4+4*i+3,1.0);
166  }
167  solver1.loadFromCoinModel(build);
168  // To make CbcBranchLink simpler assume that all variables with same i are consecutive
169 
170  double time1 = CoinCpuTime();
171  solver1.initialSolve();
172  solver1.writeMps("bad");
173  CbcModel model(solver1);
174  model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
175  // Although just one set - code as if more
176  CbcObject ** objects = new CbcObject * [1];
177  /* Format is number in sets, number in each link, first variable in matrix)
178      and then a weight for each in set to say where to branch. 
179      Finally a set number as ID.
180  */
181  double where[]={0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0};
182  objects[0]=new CbcLink(&model,11,4,4,where,0);
183  model.addObjects(1,objects);
184  delete objects[0];
185  delete [] objects;
186  // Do complete search
187 
188  model.branchAndBound();
189
190  std::cout<<"took "<<CoinCpuTime()-time1<<" seconds, "
191           <<model.getNodeCount()<<" nodes with objective "
192           <<model.getObjValue()
193           <<(!model.status() ? " Finished" : " Not finished")
194           <<std::endl;
195
196
197  if (model.getMinimizationObjValue()<1.0e50) {
198   
199    const double * solution = model.bestSolution();
200    // check correct
201    int which=-1;
202    for (int i=0;i<11;i++) {
203      for (int j=4+4*i;j<4+4*i+4;j++) {
204        double value=solution[j];
205        if (fabs(value)>1.0e-7) {
206          if (which==-1)
207            which=i;
208          else
209            assert (which==i);
210        }
211      }
212    }
213    double x=solution[0];
214    double y=solution[1];
215    double z=solution[2];
216    double w=solution[3];
217    // check z
218    assert (fabs(z-0.1*((double) which))<1.0e-7);
219    printf("Optimal solution when x is %g, y %g, z %g and w %g\n",
220           x,y,z,w);
221    printf("solution should be %g\n",10.0*x+y+z-w);
222  }
223  return 0;
224}   
Note: See TracBrowser for help on using the repository browser.