source: trunk/Cbc/examples/gear.cpp

Last change on this file was 2469, checked in by unxusr, 3 months ago

formatting

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