Last change on this file since 1574 was 1574, checked in by lou, 10 years ago

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