source: trunk/Cbc/examples/gear.cpp @ 1574

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

Change to EPL license notice.

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