Last change on this file since 2186 was 1898, checked in by stefan, 6 years ago

fixup examples

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