1 /* $Id$
2 *
3 * Name: conv-exprSinCos.cpp
4 * Author: Pietro Belotti
5 * Purpose: convexification methods for sines and cosines
6 *
7 * (C) Carnegie-Mellon University, 2006-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include <math.h>
12 #ifndef M_PI
13 # define M_PI 3.14159265358979323846
14 #endif
15 #ifndef M_PI_2
16 # define M_PI_2 1.57079632679489661923
17 #endif
18
19 #include "CouenneCutGenerator.hpp"
20
21 #include "OsiSolverInterface.hpp"
22 #include "CouenneTypes.hpp"
23 #include "CouenneProblem.hpp"
24
25 #include "CouenneExprSin.hpp"
26 #include "CouenneExprCos.hpp"
27 #include "CouenneExprAux.hpp"
28
29 namespace Couenne {
30
31 #define NEW_TRIG
32
33 #ifndef NEW_TRIG
34 /// add four cuts with slope 1 and -1
35 int addHexagon (const CouenneCutGenerator *, // cut generator that has called us
36 OsiCuts &, // cut set to be enriched
37 enum cou_trig, // sine or cosine
38 expression *, // auxiliary variable
39 expression *); // argument of cos/sin (should be a variable)
40 #endif
41
42 /// convex cuts for sine or cosine
43 int trigEnvelope (const CouenneCutGenerator *, OsiCuts &,
44 expression *, expression *, enum cou_trig);
45
46
47 /// generate convexification cut for constraint w = sin (this)
48
generateCuts(expression * w,OsiCuts & cs,const CouenneCutGenerator * cg,t_chg_bounds * chg,int wind,CouNumber lbw,CouNumber ubw)49 void exprSin::generateCuts (expression *w, //const OsiSolverInterface &si,
50 OsiCuts &cs, const CouenneCutGenerator *cg,
51 t_chg_bounds *chg, int wind,
52 CouNumber lbw, CouNumber ubw) {
53
54 // int wi = w -> Index ();
55
56 /*if (chg && !(cg -> isFirst ()) &&
57 (chg [wi].lower() == t_chg_bounds::UNCHANGED) &&
58 (chg [wi].upper() == t_chg_bounds::UNCHANGED))
59 return;*/
60
61 #ifdef NEW_TRIG
62 if (trigEnvelope (cg, cs, w, w -> Image () -> Argument (), COU_SINE) == 0)
63 #else
64 if (addHexagon (cg, cs, COU_SINE, w, w -> Image () -> Argument()) == 0)
65 #endif
66 {
67
68 }
69 }
70
71
72 /// generate convexification cut for constraint w = cos (this)
73
generateCuts(expression * w,OsiCuts & cs,const CouenneCutGenerator * cg,t_chg_bounds * chg,int wind,CouNumber lbw,CouNumber ubw)74 void exprCos::generateCuts (expression *w, //const OsiSolverInterface &si,
75 OsiCuts &cs, const CouenneCutGenerator *cg,
76 t_chg_bounds *chg, int wind,
77 CouNumber lbw, CouNumber ubw) {
78
79 // int wi = w -> Index ();
80
81 /*if (chg && !(cg -> isFirst ()) &&
82 (chg [wi].lower() == t_chg_bounds::UNCHANGED) &&
83 (chg [wi].upper() == t_chg_bounds::UNCHANGED))
84 return;*/
85
86 #ifdef NEW_TRIG
87 if (trigEnvelope (cg, cs, w, w -> Image () -> Argument (), COU_COSINE) == 0)
88 #else
89 if (addHexagon (cg, cs, COU_COSINE, w, w -> Image () -> Argument()) == 0)
90 #endif
91 {
92
93 }
94 }
95
96
97 /// restrict to quarter of the interval [0,2pi]
98 int bayEnvelope (const CouenneCutGenerator *, OsiCuts &, int, int,
99 CouNumber, CouNumber, CouNumber, bool &, bool &);
100
101
102 /// real linearization of sine/cosine
103
trigEnvelope(const CouenneCutGenerator * cg,OsiCuts & cs,expression * w,expression * arg,enum cou_trig which_trig)104 int trigEnvelope (const CouenneCutGenerator *cg, // cut generator that has called us
105 OsiCuts &cs, // cut set to be enriched
106 expression *w,
107 expression *arg,
108 enum cou_trig which_trig) {
109
110 CouNumber lb, ub;
111 arg -> getBounds (lb, ub);
112
113 // if cosine, scale variables to pretend this is a sine problem
114 CouNumber displ = (which_trig == COU_COSINE) ? M_PI_2 : 0.;
115
116 int ncuts = 0,
117 xi = arg -> Index (),
118 wi = w -> Index ();
119
120 if (fabs (ub - lb) < COUENNE_EPS) {
121
122 CouNumber x0 = 0.5 * (ub+lb), f, fp;
123
124 if (which_trig == COU_SINE) {f = sin (x0); fp = cos (x0);}
125 else {f = cos (x0); fp = -sin (x0);}
126
127 return cg -> createCut (cs, f - fp*x0, cg -> Problem () -> Var (wi) -> sign (), wi, 1., xi, -fp);
128 }
129
130 // true if, in the first call (lb), a lower/upper chord was added
131 // --> no such chord must be generated in the second call (ub)
132 bool skip_up = false,
133 skip_dn = false;
134
135 if (lb > -COUENNE_INFINITY) ncuts += bayEnvelope (cg, cs, wi, xi, lb, ub, displ, skip_up, skip_dn);
136 if (ub < COUENNE_INFINITY) ncuts += bayEnvelope (cg, cs, wi, xi, ub, lb, displ, skip_up, skip_dn);
137
138 return ncuts;
139 }
140
141
142 // __
143 // study single bay ( \__/ or / \ ) of the trigonometric function
144 //
145
bayEnvelope(const CouenneCutGenerator * cg,OsiCuts & cs,int wi,int xi,CouNumber x0,CouNumber x1,CouNumber displacement,bool & skip_up,bool & skip_dn)146 int bayEnvelope (const CouenneCutGenerator *cg, // cut generator that has called us
147 OsiCuts &cs, // cut set to be enriched
148 int wi, // dependent variable's index
149 int xi, // independent
150 CouNumber x0, // starting point
151 CouNumber x1, // other bound
152 CouNumber displacement,
153 bool &skip_up,
154 bool &skip_dn) {
155
156 enum expression::auxSign sign = cg -> Problem () -> Var (wi) -> sign ();
157
158 CouNumber tpt,
159 rx0 = modulo (x0 + displacement, 2*M_PI),
160 rx1 = rx0 + x1 - x0,
161 base = x0 - rx0,
162 sinrx0 = sin (rx0), zero;
163
164 int ncuts = 0,
165 up = (rx0 < M_PI) ? +1 : -1,
166 left = (x0 < x1) ? +1 : -1;
167
168 // starting point of the current bay
169 zero = (up>0) ? 0. : M_PI;
170
171 bool *s0, *s1;
172
173 if (up>0) {s0 = &skip_up; s1 = &skip_dn;}
174 else {s0 = &skip_dn; s1 = &skip_up;}
175
176 if (left * (modulo (rx0, M_PI) - M_PI_2) < 0) {
177
178 // after flex (i.e., at \_ or /~ ) for left bound,
179 // before flex (i.e., at _/ or ~\ ) for right bound
180
181 // out of the "belly": tangent. If on upper bay consider the lower
182 // half-plane, and viceversa --> use -up
183 if (sign != up)
184 ncuts += cg -> addTangent (cs, wi, xi, x0, sin (rx0), cos (rx0), -up);
185
186 // leftmost extreme to search for tangent point
187 CouNumber extr0 = .75 * M_PI * (left+1) - M_PI_2 * up;
188
189 // in:
190 if ((left * (rx1 - M_PI * ((left - up) / 2 + 1)) <= 0) || // if rx1 in same "belly", or
191 (left * (rx1 - (tpt = trigNewton
192 (rx0, extr0, extr0 + M_PI_2))) <= 0)) { // before closest leaning point
193 if (!*s1 && (sign != -up)) // -> chord, if not already added in previous call
194 *s1 = ((ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1, sin (rx1), up)) > 0);
195 } else
196 if (sign != -up)
197 ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base+tpt, sin (tpt), up);
198 } else {
199
200 // after stationary point (i.e., _/ or ~\ ) for left bound,
201 // before stationary point (i.e., /~ or \_ ) for right bound
202
203 // if (left * (rx1 - left * (zero + 5*M_PI_2)) < 0) {
204 if (left * (rx1 - (4*left - up + 2) * M_PI_2) < 0) {
205 CouNumber cosrx0 = cos (rx0);
206 if (up * (sin (rx1) - sinrx0 - cosrx0 * (rx1-rx0)) < 0) {
207 // (b,sinb) below tangent --> tangent
208 if (sign != up)
209 ncuts += cg -> addTangent (cs, wi, xi, x0, sinrx0, cosrx0, -up);
210 } else { // up: either chord or leaning plane
211 CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
212 tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
213 if (left * (rx1 - tpt) < 0) {
214 if (!*s0 && (sign != up) )
215 *s0 = ((ncuts += cg->addSegment (cs, wi, xi, x0, sin(rx0), x1, sin(rx1), -up)) > 0);
216 } else
217 if (sign != up)
218 ncuts += cg->addSegment (cs, wi, xi, x0, sin(rx0), base+tpt, sin(tpt), -up);
219 }
220 } else {
221 CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
222 tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
223 if (sign != up)
224 ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), -up);
225 }
226
227 // down: other chord or leaning plane
228 if ((left * (rx1 - (zero + M_PI)) < 0) ||
229 (left * (rx1 - (tpt = trigNewton (rx0, (2 + left - up) * M_PI_2,
230 (2 + 2*left - up) * M_PI_2))) < 0)) {
231 if (!*s1 && (sign != -up))
232 *s1 = ((ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1, sin (rx1), up)) > 0);
233 } else
234 if (sign != -up)
235 ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), up);
236 }
237
238 return ncuts;
239 }
240
241
242 #ifndef NEW_TRIG
243
244
245 /// add lateral edges of the hexagon providing
246
addHexagon(const CouenneCutGenerator * cg,OsiCuts & cs,enum cou_trig tt,expression * aux,expression * arg)247 int addHexagon (const CouenneCutGenerator *cg, // cut generator that has called us
248 OsiCuts &cs, // cut set to be enriched
249 enum cou_trig tt, // sine or cosine
250 expression *aux, // auxiliary variable
251 expression *arg) { // argument of cos/sin (should be a variable)
252
253 // retrieve argument bounds
254 CouNumber lb, ub;
255 arg -> getBounds (lb, ub);
256
257 int ncuts = 0,
258 x_ind = arg -> Index (),
259 w_ind = aux -> Index ();
260
261 enum auxSign sign = cg -> Problem () -> Var (w_ind) -> sign ();
262
263 if (fabs (ub - lb) < COUENNE_EPS) {
264
265 CouNumber x0 = 0.5 * (ub+lb), f, fp;
266
267 if (tt == COU_SINE) {f = sin (x0); fp = cos (x0);}
268 else {f = cos (x0); fp = -sin (x0);}
269
270 return cg -> createCut (cs, f - fp*x0, sign, w_ind, 1., x_ind, -fp);
271 }
272
273 // add / \ envelope
274 // \ /
275
276 // left
277 if (lb > -COUENNE_INFINITY) { // if not unbounded
278 if (tt == COU_SINE) {
279 if (sign != expression::AUX_GEQ) ncuts += cg -> createCut (cs, sin (lb) - lb, -1, w_ind, 1., x_ind, -1.); // up: w-x <= f lb - lb
280 if (sign != expression::AUX_LEQ) ncuts += cg -> createCut (cs, sin (lb) + lb, +1, w_ind, 1., x_ind, 1.); // dn: w+x >= f lb + lb
281 }
282 else {
283 if (sign != expression::AUX_GEQ) ncuts += cg -> createCut (cs, cos (lb) - lb, -1, w_ind, 1., x_ind, -1.); // up: w-x <= f lb - lb
284 if (sign != expression::AUX_LEQ) ncuts += cg -> createCut (cs, cos (lb) + lb, +1, w_ind, 1., x_ind, 1.); // dn: w+x >= f lb + lb
285 }
286 }
287
288 // right
289 if (ub < COUENNE_INFINITY) { // if not unbounded
290 if (tt == COU_SINE) {
291 if (sign != expression::AUX_LEQ) ncuts += cg -> createCut (cs, sin (ub) - ub, +1, w_ind, 1., x_ind, -1.); // dn: w-x >= f ub - ub
292 if (sign != expression::AUX_GEQ) ncuts += cg -> createCut (cs, sin (ub) + ub, -1, w_ind, 1., x_ind, 1.); // up: w+x <= f ub + ub
293 }
294 else {
295 if (sign != expression::AUX_LEQ) ncuts += cg -> createCut (cs, cos (ub) - ub, +1, w_ind, 1., x_ind, -1.); // dn: w-x >= f ub - ub
296 if (sign != expression::AUX_GEQ) ncuts += cg -> createCut (cs, cos (ub) + ub, -1, w_ind, 1., x_ind, 1.); // up: w+x <= f ub + ub
297 }
298 }
299
300 return ncuts;
301 }
302
303 #endif
304
305 }
306