source: trunk/Cbc/examples/link.cpp @ 2469

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

formatting

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