source: branches/Couenne/Couenne/src/convex/operators/conv-exprSinCos.cpp @ 534

Last change on this file since 534 was 534, checked in by pbelotti, 13 years ago

moved include files to make them doxygenable. Introduced three-way branching, with fixed intervals for now. Added check for small bound interval within all generateCuts()

File size: 7.1 KB
Line 
1/*
2 * Name:    conv-exprSinCos.cpp
3 * Author:  Pietro Belotti
4 * Purpose: convexification methods for sines and cosines
5 *
6 * This file is licensed under the Common Public License (CPL)
7 */
8
9#include <math.h>
10
11#include <OsiSolverInterface.hpp>
12#include <CouenneTypes.h>
13#include <CouenneCutGenerator.h>
14#include <exprSin.h>
15#include <exprCos.h>
16#include <exprAux.h>
17
18//#define NEW_TRIG
19
20/// convex cuts for sine or cosine
21void trigEnvelope (const CouenneCutGenerator *, OsiCuts &,
22                   exprAux *, expression *, enum cou_trig);
23
24
25// generate convexification cut for constraint w = sin (this)
26
27void exprSin::generateCuts (exprAux *w, const OsiSolverInterface &si, 
28                            OsiCuts &cs, const CouenneCutGenerator *cg) {
29
30#ifdef NEW_TRIG
31  trigEnvelope (cg, cs, w, w -> Image () -> Argument (), COU_SINE);
32#else
33  addHexagon (cg, cs, sin, w, w -> Image () -> Argument());
34#endif
35}
36
37
38// generate convexification cut for constraint w = cos (this)
39
40void exprCos::generateCuts (exprAux *w, const OsiSolverInterface &si, 
41                            OsiCuts &cs, const CouenneCutGenerator *cg) {
42
43#ifdef NEW_TRIG
44  trigEnvelope (cg, cs, w, w -> Image () -> Argument (), COU_COSINE);
45#else
46  addHexagon (cg, cs, cos, w, w -> Image () -> Argument());
47#endif
48}
49
50
51// add lateral edges of the hexagon providing
52
53void addHexagon (const CouenneCutGenerator *cg, // cut generator that has called us
54                 OsiCuts &cs,       // cut set to be enriched
55                 unary_function f,  // sine or cosine
56                 exprAux *aux,      // auxiliary variable
57                 expression *arg) { // argument of cos/sin (should be a variable)
58
59  expression *lbe, *ube;
60  arg -> getBounds (lbe, ube);
61
62  CouNumber lb = (*lbe) (), 
63            ub = (*ube) ();
64
65  int x_ind = arg -> Index ();
66  int w_ind = aux -> Index ();
67  /*
68  if (fabs (ub - lb) < COUENNE_EPS) {
69
70    CouNumber x0 = 0.5 * (ub+lb), f, fp;
71
72    if (which_trig == COU_SINE) {f = sin (x0); fp =  cos (x0);}
73    else                        {f = cos (x0); fp = -sin (x0);}
74
75    cg -> createCut (cs, f - fp*x0, 0, w_ind, 1., x_ind, -fp);
76    return;
77  }
78  */
79  // add the lower envelope
80  cg -> createCut (cs, f (lb) - lb, -1, w_ind, 1., x_ind, -1.); // left:  w - x <= f lb - lb
81  cg -> createCut (cs, f (ub) + ub, -1, w_ind, 1., x_ind,  1.); // right: w + x <= f ub + ub
82
83  // add the upper envelope
84  cg -> createCut (cs, f (ub) - ub, +1, w_ind, 1., x_ind, -1.); // right: w - x >= cos ub - ub
85  cg -> createCut (cs, f (lb) + lb, +1, w_ind, 1., x_ind,  1.); // left:  w + x >= cos lb + lb
86
87  delete lbe;
88  delete ube;
89}
90
91
92/// normalize angle within [0,b] (typically, pi or 2pi)
93inline CouNumber modulo (register CouNumber a, register CouNumber b)
94  {return a - b * floor (a/b);}
95
96
97/// restrict to quarter of the interval [0,2pi]
98int bayEnvelope (const CouenneCutGenerator *, OsiCuts &, int, int, 
99                 CouNumber, CouNumber, CouNumber, bool &, bool &);
100
101
102/// real linearization of sine/cosine
103
104void trigEnvelope (const CouenneCutGenerator *cg, // cut generator that has called us
105                   OsiCuts &cs,                   // cut set to be enriched
106                   exprAux *w,
107                   expression *arg,
108                   enum cou_trig which_trig) {
109
110  expression *lbe, *ube;
111
112  arg -> getBounds (lbe, ube);
113
114  CouNumber lb = (*lbe) (), 
115            ub = (*ube) (),
116            // if cosine, scale variables to pretend this is a sine problem
117            displacement = (which_trig == COU_COSINE) ? M_PI_2 : 0.;
118
119  delete lbe;
120  delete ube;
121
122  int xi = arg -> Index (),
123      wi = w   -> Index ();
124
125  if (fabs (ub - lb) < COUENNE_EPS) {
126
127    CouNumber x0 = 0.5 * (ub+lb), f, fp;
128
129    if (which_trig == COU_SINE) {f = sin (x0); fp =  cos (x0);}
130    else                        {f = cos (x0); fp = -sin (x0);}
131
132    cg -> createCut (cs, f - fp*x0, 0, wi, 1., xi, -fp);
133    return;
134  }
135
136  // true if, in the first call (lb), a lower/upper chord was added
137  // --> no such chord must be generated in the second call (ub)
138  bool skip_up = false, 
139       skip_dn = false;
140
141  bayEnvelope (cg, cs, wi, xi, lb, ub, displacement, skip_up, skip_dn); // lb
142  bayEnvelope (cg, cs, wi, xi, ub, lb, displacement, skip_up, skip_dn); // ub
143}
144
145
146// study single bay ( \__/ or /~~\ ) of the trigonometric function
147
148int bayEnvelope (const CouenneCutGenerator *cg, // cut generator that has called us
149                 OsiCuts &cs,                   // cut set to be enriched
150                 int wi,
151                 int xi,
152                 CouNumber x0,           // starting point
153                 CouNumber x1,           // other bound
154                 CouNumber displacement, 
155                 bool &skip_up, 
156                 bool &skip_dn) {
157
158  CouNumber tpt,
159    rx0  = modulo (x0 + displacement, 2*M_PI),
160    rx1  = rx0 + x1 - x0,
161    base = x0 + displacement - rx0,
162    sinrx0 = sin (rx0), zero;
163
164  int
165    up   = (modulo (rx0, 2*M_PI) < M_PI) ? +1 : -1,
166    left = (x0 < x1)                     ? +1 : -1;
167
168  zero = (up>0) ? 0. : M_PI;
169
170  bool *s0, *s1;
171
172  if (up>0) {s0 = &skip_up; s1 = &skip_dn;}
173  else      {s0 = &skip_dn; s1 = &skip_up;}
174
175  if (left * (modulo (rx0, M_PI) - M_PI_2) < 0) { 
176
177    // after  flex (i.e., at \_ or /~ ) for left  bound,
178    // before flex (i.e., at _/ or ~\ ) for right bound
179
180    // out of the "belly": tangent. If on upper bay we consider the
181    // lower half-plane, and viceversa --> use -up
182    cg -> addTangent (cs, wi, xi, x0, sin (rx0), cos (rx0), -up);
183
184    // leftmost extreme to search for tangent point
185    CouNumber extr0 = .75 * M_PI * (left+1) - M_PI_2 * up; 
186
187    // in:
188    if ((left * (rx1 - M_PI * ((left - up) / 2 + 1)) <= 0) ||   // if rx1 in same "belly", or
189        (left * (rx1 - (tpt = trigNewton
190                        (rx0, extr0, extr0 + M_PI_2))) <= 0)) { // before closest leaning point
191      if (!*s1) // -> chord, if not already added in previous call
192        *s1 = (cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1,         sin (rx1), up) > 0);
193    } else     cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), up);
194  }
195  else {
196
197    // after  stationary point (i.e., _/ or ~\ ) for left  bound,
198    // before stationary point (i.e., /~ or \_ ) for right bound
199 
200    //    if (left * (rx1 - left * (zero + 5*M_PI_2)) < 0) {
201    if (left * (rx1 - (4*left - up + 2) * M_PI_2) < 0) {
202      CouNumber cosrx0 = cos (rx0);
203      if (up * (sin (rx1) - sinrx0 - cosrx0 * (rx1-rx0)) < 0) 
204        // (b,sinb) below tangent --> tangent
205        cg -> addTangent (cs, wi, xi, x0, sinrx0, cosrx0, -up);
206      else {    // up: either chord or leaning plane
207        CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
208        tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
209        if (left * (rx1 - tpt) < 0) {
210          if (!*s0)
211            *s0 = cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1,         sin (rx1), -up) > 0;
212        } else    cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), -up);
213      }
214    } else {
215      CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
216      tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
217      cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), -up);
218    }
219
220    // down: other chord or leaning plane
221    if ((left * (rx1 - (zero + M_PI)) < 0) || 
222        (left * (rx1 - (tpt = trigNewton (rx0, 
223                                          (2 +   left - up) * M_PI_2, 
224                                          (2 + 2*left - up) * M_PI_2))) < 0)) {
225      if (!*s1) 
226        *s1 = (cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1,         sin (rx1), up) > 0);
227    } else     cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), up);
228  }
229
230  return 0;
231}
Note: See TracBrowser for help on using the repository browser.