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