source: branches/devel/Cbc/examples/gear.cpp @ 425

Last change on this file since 425 was 128, checked in by forrest, 15 years ago

add reference

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.2 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  abs  ( 1.0/6.931 - x1*x4/x2*x3)
31
32where the variables are integral between 12 and 60.
33See E.Sangren, "Nonlinear Integer and Discrete Programming in
34Mechanical Design Optimization". Trans. ASME, J. Mech Design 112, 223-229, 1990
35
36One could try to use logarithms to make the problem separable but that leads to a
37weak formulation.  Instaed we are going to use linked
38special ordered sets.  The generalization with column generation can be even more powerful
39but is not yet in CBC.
40
41The idea is simple:
42
43A linear variable is a convex combination of its lower bound and upper bound!
44If x must lie between 12 and 60 then we can substitute for x  as x == 12.0*xl + 60.0*xu where
45xl + xu == 1.0.  At first this looks cumbersome but if we have xl12, xl13, ... xl60 and corresponding
46xu and yl and yu then we can write:
47
48x == sum 12.0*xl[i] + 60.0* xu[i] where sum xl[i] + xu[i] == 1.0
49and
50x*y == 12.0*12.0*xl12 + 12.0*60.0*xu12 + 13.0*12.0*xl13 + 13.0*60.0*x13 ....
51                 + 12.0*60*.0xl60 + 60.0*60.0*xu60
52
53And now x*y is correct if x is integer and xl[i], xu[i] are only nonzero for one i.
54Note that this would have worked just as easily for y**2 or any clean function of y.
55
56So this is just like a special ordered set of type 1 but on two sets simultaneously.
57The idea is even more powerful if we want other functions on y as we can branch on all
58sets 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  /*
90    We are going to treat x1 and x2 as integer and x3 and x4 as a set.
91    We define two new variables y1 == x1*x4 and y2 == x2*x3.
92    We define a variable z == x1*x4/x2*x3 so y2*z == y1
93    (we will treat y2 as a set)
94    Then we have objective - minimize w1 + w2 where
95    w1 - w2 = 1.0/6.931 - z
96
97    The model would be a lot smaller if we had column generation.
98  */
99  // Create model
100  CoinModel build;
101  // Keep values of all variables for reporting purposes even if not necessary
102  /*
103    z is first, then x then y1,y2 then w1,w2
104    then y1 stuff, y2 stuff and finally y2 -> z stuff.
105    For rows same but 2 per y then rest of z stuff
106  */
107  int loInt=12;
108  int hiInt=60;
109  int ybaseA=5, ybaseB=9, ylen=hiInt-loInt+1;
110  int base = ybaseB+2*2*ylen;
111  int yylen = hiInt*hiInt-loInt*loInt+1;
112  int zbase = 10;
113  int i;
114  // Do single variables
115  double value[] ={1.0,1.0};
116  int row[2];
117  /* z - obviously we can't choose bounds too tight but we need bounds
118     so choose 20% off as obviously feasible.
119     fastest way to solve would be too run for a few seconds to get
120     tighter bounds then re-formulate and solve.  */
121  double loose=0.2;
122  double loZ = (1-loose)*(1.0/6.931), hiZ = (1+loose)*(1.0/6.931);
123  row[0]=0; // for reporting
124  row[1]=zbase+1; // for real use
125  build.addColumn(2,row,value,loZ, hiZ, 0.0);
126  // x
127  for (i=0;i<4;i++) {
128    row[0]=i+1;
129    build.addColumn(1,row,value,loInt, hiInt,0.0);
130    // we don't need to say x2, x3 integer but won't hurt
131    build.setInteger(i+1);
132  }
133  // y
134  for (i=0;i<2;i++) {
135    // y from x*x, and convexity
136    row[0]=ybaseA+2*i;
137    if (i==0)
138      row[1]=zbase+2; // yb*z == ya
139    else
140      row[1]=zbase-1; // to feed into z
141    build.addColumn(2,row,value,loInt*loInt, hiInt*hiInt,0.0);
142    // we don't need to say integer but won't hurt
143    build.setInteger(ybaseA+i);
144  }
145  // skip z convexity put w in final equation
146  row[0]=zbase+1;
147  build.addColumn(1,row,value,0.0,1.0,1.0);
148  value[0]=-1.0;
149  build.addColumn(1,row,value,0.0,1.0,1.0);
150  // Do columns so we know where each is
151  for (i=ybaseB;i<base+(2*yylen);i++)
152    build.setColumnBounds(i,0.0,1.0);
153  // Now do rows
154  // z definition
155  build.setRowBounds(0,0.0,0.0);
156  for (i=0;i<yylen;i++) {
157    // l
158    build.setElement(0,base+2*i,-loZ);
159    // u
160    build.setElement(0,base+2*i+1,-hiZ);
161  }
162  // x
163  for (i=0;i<2;i++) {
164    int iVarRow = 1+i;
165    int iSetRow = 4-i; // as it is x1*x4 and x2*x3
166    build.setRowBounds(iVarRow,0.0,0.0);
167    build.setRowBounds(iSetRow,0.0,0.0);
168    int j;
169    int base2 = ybaseB + 2*ylen*i;
170    for (j=0;j<ylen;j++) {
171      // l
172      build.setElement(iVarRow,base2+2*j,-loInt);
173      build.setElement(iSetRow,base2+2*j,-loInt-j);
174      // u
175      build.setElement(iVarRow,base2+2*j+1,-hiInt);
176      build.setElement(iSetRow,base2+2*j+1,-loInt-j);
177    }
178  }
179  // y
180  for (i=0;i<2;i++) {
181    int iRow = 5+2*i;
182    int iConvex = iRow+1;
183    build.setRowBounds(iRow,0.0,0.0);
184    build.setRowBounds(iConvex,1.0,1.0);
185    int j;
186    int base2 = ybaseB + 2*ylen*i;
187    for (j=0;j<ylen;j++) {
188      // l
189      build.setElement(iRow,base2+2*j,-loInt*(j+loInt));
190      build.setElement(iConvex,base2+2*j,1.0);
191      // u
192      build.setElement(iRow,base2+2*j+1,-hiInt*(j+loInt));
193      build.setElement(iConvex,base2+2*j+1,1.0);
194    }
195  }
196  // row that feeds into z and convexity
197  build.setRowBounds(zbase-1,0.0,0.0);
198  build.setRowBounds(zbase,1.0,1.0);
199  for (i=0;i<yylen;i++) {
200    // l
201    build.setElement(zbase-1,base+2*i,-(i+loInt*loInt));
202    build.setElement(zbase,base+2*i,1.0);
203    // u
204    build.setElement(zbase-1,base+2*i+1,-(i+loInt*loInt));
205    build.setElement(zbase,base+2*i+1,1.0);
206  }
207  // and real equation rhs
208  build.setRowBounds(zbase+1,1.0/6.931,1.0/6.931);
209  // z*y
210  build.setRowBounds(zbase+2,0.0,0.0);
211  for (i=0;i<yylen;i++) {
212    // l
213    build.setElement(zbase+2,base+2*i,-(i+loInt*loInt)*loZ);
214    // u
215    build.setElement(zbase+2,base+2*i+1,-(i+loInt*loInt)*hiZ);
216  }
217  // And finally two more rows to break symmetry
218  build.setRowBounds(zbase+3,-COIN_DBL_MAX,0.0);
219  build.setElement(zbase+3,1,1.0);
220  build.setElement(zbase+3,4,-1.0);
221  build.setRowBounds(zbase+4,-COIN_DBL_MAX,0.0);
222  build.setElement(zbase+4,2,1.0);
223  build.setElement(zbase+4,3,-1.0);
224  solver1.loadFromCoinModel(build);
225  // To make CbcBranchLink simpler assume that all variables with same i are consecutive
226 
227  double time1 = CoinCpuTime();
228  solver1.initialSolve();
229  solver1.writeMps("bad");
230  CbcModel model(solver1);
231  model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
232  model.solver()->setHintParam(OsiDoScale,false,OsiHintTry);
233
234  CbcObject ** objects = new CbcObject * [3];
235  /* Format is number in sets, number in each link, first variable in matrix)
236      and then a weight for each in set to say where to branch. 
237      In this case use NULL to say 0,1,2 ...
238      Finally a set number as ID.
239  */
240  objects[0]=new CbcLink(&model,ylen,2,ybaseB,NULL,0);
241  objects[0]->setPriority(10);
242  objects[1]=new CbcLink(&model,ylen,2,ybaseB+2*ylen,NULL,0);
243  objects[1]->setPriority(20);
244  objects[2]=new CbcLink(&model,yylen,2,base,NULL,0);
245  objects[2]->setPriority(1);
246  model.addObjects(3,objects);
247  for (i=0;i<3;i++)
248    delete objects[i];
249  delete [] objects;
250  model.messageHandler()->setLogLevel(1);
251  // Do complete search
252 
253  model.setDblParam(CbcModel::CbcMaximumSeconds,1200.0);
254  model.setDblParam(CbcModel::CbcCutoffIncrement,1.0e-8);
255  model.branchAndBound();
256
257  std::cout<<"took "<<CoinCpuTime()-time1<<" seconds, "
258           <<model.getNodeCount()<<" nodes with objective "
259           <<model.getObjValue()
260           <<(!model.status() ? " Finished" : " Not finished")
261           <<std::endl;
262
263
264  if (model.getMinimizationObjValue()<1.0e50) {
265   
266    const double * solution = model.bestSolution();
267    int numberColumns = model.solver()->getNumCols();
268    double x1=solution[1];
269    double x2=solution[2];
270    double x3=solution[3];
271    double x4=solution[4];
272    printf("Optimal solution %g %g %g %g\n",x1,x2,x3,x4);
273    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
274      double value=solution[iColumn];
275      if (fabs(value)>1.0e-7) 
276        std::cout<<iColumn<<" "<<value<<std::endl;
277    }
278  }
279  return 0;
280}   
Note: See TracBrowser for help on using the repository browser.