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

Last change on this file since 1898 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: 9.2 KB
Line 
1// $Id: gear.cpp 1898 2013-04-09 18:06:04Z stefan $
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#include <cassert>
7#include <iomanip>
8
9#include "CoinPragma.hpp"
10
11// For Branch and bound
12#include "OsiSolverInterface.hpp"
13#include "CbcModel.hpp"
14#include "CoinModel.hpp"
15// For Linked Ordered Sets
16#include "CbcBranchLink.hpp"
17#include "OsiClpSolverInterface.hpp"
18
19#include "CoinTime.hpp"
20
21
22/************************************************************************
23
24This shows how we can define a new branching method to solve problems with
25nonlinearities and discontinuities.
26
27We are going to solve the problem
28
29minimize  abs  ( 1.0/6.931 - x1*x4/x2*x3)
30
31where the variables are integral between 12 and 60.
32See E.Sangren, "Nonlinear Integer and Discrete Programming in
33Mechanical Design Optimization". Trans. ASME, J. Mech Design 112, 223-229, 1990
34
35One could try to use logarithms to make the problem separable but that leads to a
36weak formulation.  Instaed we are going to use linked
37special ordered sets.  The generalization with column generation can be even more powerful
38but is not yet in CBC.
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 12 and 60 then we can substitute for x  as x == 12.0*xl + 60.0*xu where
44xl + xu == 1.0.  At first this looks cumbersome but if we have xl12, xl13, ... xl60 and corresponding
45xu and yl and yu then we can write:
46
47x == sum 12.0*xl[i] + 60.0* xu[i] where sum xl[i] + xu[i] == 1.0
48and
49x*y == 12.0*12.0*xl12 + 12.0*60.0*xu12 + 13.0*12.0*xl13 + 13.0*60.0*x13 ....
50                 + 12.0*60*.0xl60 + 60.0*60.0*xu60
51
52And now x*y is correct if x is integer and xl[i], xu[i] are only nonzero for one i.
53Note that this would have worked just as easily for y**2 or any clean function of y.
54
55So this is just like a special ordered set of type 1 but on two sets simultaneously.
56The idea is even more powerful if we want other functions on y as we can branch on all
57sets simultaneously.
58Also note that convexity requirements for any non-linear functions are not needed.
59
60So we need a new branching method to do that - see CbcBranchLink.?pp
61
62We are going to need a CbcBranchLink method to see whether we are satisfied etc and also to
63create another branching object which knows how to fix variables.  We might be able to use an
64existing method for the latter but let us create two methods CbcLink and
65CbcLinkBranchingObject.
66
67For CbcLink we will need the following methods:
68Constructot/Destructor
69infeasibility - returns 0.0 if feasible otherwise some measure of infeasibility
70feasibleRegion - sets bounds to contain current solution
71createBranch - creates a CbcLinkBranchingObject
72
73For CbcLinkBranchingObject we need:
74Constructor/Destructor
75branch - does actual fixing
76print - optional for debug purposes.
77
78The easiest way to do this is to cut and paste from CbcBranchActual to get current
79SOS stuff and then modify that.
80
81************************************************************************/
82
83int main (int argc, const char *argv[])
84{
85
86  OsiClpSolverInterface solver1;
87
88  /*
89    We are going to treat x1 and x2 as integer and x3 and x4 as a set.
90    We define two new variables y1 == x1*x4 and y2 == x2*x3.
91    We define a variable z == x1*x4/x2*x3 so y2*z == y1
92    (we will treat y2 as a set)
93    Then we have objective - minimize w1 + w2 where
94    w1 - w2 = 1.0/6.931 - z
95
96    The model would be a lot smaller if we had column generation.
97  */
98  // Create model
99  CoinModel build;
100  // Keep values of all variables for reporting purposes even if not necessary
101  /*
102    z is first, then x then y1,y2 then w1,w2
103    then y1 stuff, y2 stuff and finally y2 -> z stuff.
104    For rows same but 2 per y then rest of z stuff
105  */
106  int loInt=12;
107  int hiInt=60;
108  int ybaseA=5, ybaseB=9, ylen=hiInt-loInt+1;
109  int base = ybaseB+2*2*ylen;
110  int yylen = hiInt*hiInt-loInt*loInt+1;
111  int zbase = 10;
112  int i;
113  // Do single variables
114  double value[] ={1.0,1.0};
115  int row[2];
116  /* z - obviously we can't choose bounds too tight but we need bounds
117     so choose 20% off as obviously feasible.
118     fastest way to solve would be too run for a few seconds to get
119     tighter bounds then re-formulate and solve.  */
120  double loose=0.2;
121  double loZ = (1-loose)*(1.0/6.931), hiZ = (1+loose)*(1.0/6.931);
122  row[0]=0; // for reporting
123  row[1]=zbase+1; // for real use
124  build.addColumn(2,row,value,loZ, hiZ, 0.0);
125  // x
126  for (i=0;i<4;i++) {
127    row[0]=i+1;
128    build.addColumn(1,row,value,loInt, hiInt,0.0);
129    // we don't need to say x2, x3 integer but won't hurt
130    build.setInteger(i+1);
131  }
132  // y
133  for (i=0;i<2;i++) {
134    // y from x*x, and convexity
135    row[0]=ybaseA+2*i;
136    if (i==0)
137      row[1]=zbase+2; // yb*z == ya
138    else
139      row[1]=zbase-1; // to feed into z
140    build.addColumn(2,row,value,loInt*loInt, hiInt*hiInt,0.0);
141    // we don't need to say integer but won't hurt
142    build.setInteger(ybaseA+i);
143  }
144  // skip z convexity put w in final equation
145  row[0]=zbase+1;
146  build.addColumn(1,row,value,0.0,1.0,1.0);
147  value[0]=-1.0;
148  build.addColumn(1,row,value,0.0,1.0,1.0);
149  // Do columns so we know where each is
150  for (i=ybaseB;i<base+(2*yylen);i++)
151    build.setColumnBounds(i,0.0,1.0);
152  // Now do rows
153  // z definition
154  build.setRowBounds(0,0.0,0.0);
155  for (i=0;i<yylen;i++) {
156    // l
157    build.setElement(0,base+2*i,-loZ);
158    // u
159    build.setElement(0,base+2*i+1,-hiZ);
160  }
161  // x
162  for (i=0;i<2;i++) {
163    int iVarRow = 1+i;
164    int iSetRow = 4-i; // as it is x1*x4 and x2*x3
165    build.setRowBounds(iVarRow,0.0,0.0);
166    build.setRowBounds(iSetRow,0.0,0.0);
167    int j;
168    int base2 = ybaseB + 2*ylen*i;
169    for (j=0;j<ylen;j++) {
170      // l
171      build.setElement(iVarRow,base2+2*j,-loInt);
172      build.setElement(iSetRow,base2+2*j,-loInt-j);
173      // u
174      build.setElement(iVarRow,base2+2*j+1,-hiInt);
175      build.setElement(iSetRow,base2+2*j+1,-loInt-j);
176    }
177  }
178  // y
179  for (i=0;i<2;i++) {
180    int iRow = 5+2*i;
181    int iConvex = iRow+1;
182    build.setRowBounds(iRow,0.0,0.0);
183    build.setRowBounds(iConvex,1.0,1.0);
184    int j;
185    int base2 = ybaseB + 2*ylen*i;
186    for (j=0;j<ylen;j++) {
187      // l
188      build.setElement(iRow,base2+2*j,-loInt*(j+loInt));
189      build.setElement(iConvex,base2+2*j,1.0);
190      // u
191      build.setElement(iRow,base2+2*j+1,-hiInt*(j+loInt));
192      build.setElement(iConvex,base2+2*j+1,1.0);
193    }
194  }
195  // row that feeds into z and convexity
196  build.setRowBounds(zbase-1,0.0,0.0);
197  build.setRowBounds(zbase,1.0,1.0);
198  for (i=0;i<yylen;i++) {
199    // l
200    build.setElement(zbase-1,base+2*i,-(i+loInt*loInt));
201    build.setElement(zbase,base+2*i,1.0);
202    // u
203    build.setElement(zbase-1,base+2*i+1,-(i+loInt*loInt));
204    build.setElement(zbase,base+2*i+1,1.0);
205  }
206  // and real equation rhs
207  build.setRowBounds(zbase+1,1.0/6.931,1.0/6.931);
208  // z*y
209  build.setRowBounds(zbase+2,0.0,0.0);
210  for (i=0;i<yylen;i++) {
211    // l
212    build.setElement(zbase+2,base+2*i,-(i+loInt*loInt)*loZ);
213    // u
214    build.setElement(zbase+2,base+2*i+1,-(i+loInt*loInt)*hiZ);
215  }
216  // And finally two more rows to break symmetry
217  build.setRowBounds(zbase+3,-COIN_DBL_MAX,0.0);
218  build.setElement(zbase+3,1,1.0);
219  build.setElement(zbase+3,4,-1.0);
220  build.setRowBounds(zbase+4,-COIN_DBL_MAX,0.0);
221  build.setElement(zbase+4,2,1.0);
222  build.setElement(zbase+4,3,-1.0);
223  solver1.loadFromCoinModel(build);
224  // To make CbcBranchLink simpler assume that all variables with same i are consecutive
225 
226  double time1 = CoinCpuTime();
227  solver1.initialSolve();
228  solver1.writeMps("bad");
229  CbcModel model(solver1);
230  model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
231  model.solver()->setHintParam(OsiDoScale,false,OsiHintTry);
232
233  CbcObject ** objects = new CbcObject * [3];
234  /* Format is number in sets, number in each link, first variable in matrix)
235      and then a weight for each in set to say where to branch. 
236      In this case use NULL to say 0,1,2 ...
237      Finally a set number as ID.
238  */
239  objects[0]=new CbcLink(&model,ylen,2,ybaseB,NULL,0);
240  objects[0]->setPriority(10);
241  objects[1]=new CbcLink(&model,ylen,2,ybaseB+2*ylen,NULL,0);
242  objects[1]->setPriority(20);
243  objects[2]=new CbcLink(&model,yylen,2,base,NULL,0);
244  objects[2]->setPriority(1);
245  model.addObjects(3,objects);
246  for (i=0;i<3;i++)
247    delete objects[i];
248  delete [] objects;
249  model.messageHandler()->setLogLevel(1);
250  // Do complete search
251 
252  model.setDblParam(CbcModel::CbcMaximumSeconds,1200.0);
253  model.setDblParam(CbcModel::CbcCutoffIncrement,1.0e-8);
254  model.branchAndBound();
255
256  std::cout<<"took "<<CoinCpuTime()-time1<<" seconds, "
257           <<model.getNodeCount()<<" nodes with objective "
258           <<model.getObjValue()
259           <<(!model.status() ? " Finished" : " Not finished")
260           <<std::endl;
261
262
263  if (model.getMinimizationObjValue()<1.0e50) {
264   
265    const double * solution = model.bestSolution();
266    int numberColumns = model.solver()->getNumCols();
267    double x1=solution[1];
268    double x2=solution[2];
269    double x3=solution[3];
270    double x4=solution[4];
271    printf("Optimal solution %g %g %g %g\n",x1,x2,x3,x4);
272    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
273      double value=solution[iColumn];
274      if (fabs(value)>1.0e-7) 
275        std::cout<<iColumn<<" "<<value<<std::endl;
276    }
277  }
278  return 0;
279}   
Note: See TracBrowser for help on using the repository browser.