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

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

introducing better bounds for sin/cos. Tool for drawing cuts

File size: 6.5 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  // add the lower envelope
69  cg -> createCut (cs, f (lb) - lb, -1, w_ind, 1., x_ind, -1.); // left:  w - x <= f lb - lb
70  cg -> createCut (cs, f (ub) + ub, -1, w_ind, 1., x_ind,  1.); // right: w + x <= f ub + ub
71
72  // add the upper envelope
73  cg -> createCut (cs, f (ub) - ub, +1, w_ind, 1., x_ind, -1.); // right: w - x >= cos ub - ub
74  cg -> createCut (cs, f (lb) + lb, +1, w_ind, 1., x_ind,  1.); // left:  w + x >= cos lb + lb
75
76  delete lbe;
77  delete ube;
78}
79
80
81/// normalize angle within [0,b] (typically, pi or 2pi)
82inline CouNumber modulo (register CouNumber a, register CouNumber b)
83  {return a - b * floor (a/b);}
84
85
86/// restrict to quarter of the interval [0,2pi]
87int bayEnvelope (const CouenneCutGenerator *, OsiCuts &, int, int, 
88                 CouNumber, CouNumber, CouNumber, bool &, bool &);
89
90
91/// real linearization of sine/cosine
92
93void trigEnvelope (const CouenneCutGenerator *cg, // cut generator that has called us
94                   OsiCuts &cs,                   // cut set to be enriched
95                   exprAux *w,
96                   expression *arg,
97                   enum cou_trig which_trig) {
98
99  expression *lbe, *ube;
100
101  arg -> getBounds (lbe, ube);
102
103  CouNumber lb = (*lbe) (), 
104            ub = (*ube) (),
105            // if cosine, scale variables to pretend this is a sine problem
106            displacement = (which_trig == COU_COSINE) ? M_PI_2 : 0.;
107
108  delete lbe;
109  delete ube;
110
111  int xi = arg -> Index (),
112      wi = w   -> Index ();
113
114  // true if, in the first call (lb), a lower/upper chord was added
115  // --> no such chord must be generated in the second call (ub)
116  bool skip_up = false, 
117       skip_dn = false;
118
119  bayEnvelope (cg, cs, wi, xi, lb, ub, displacement, skip_up, skip_dn); // lb
120  bayEnvelope (cg, cs, wi, xi, ub, lb, displacement, skip_up, skip_dn); // ub
121}
122
123
124// study single bay ( \__/ or /~~\ ) of the trigonometric function
125
126int bayEnvelope (const CouenneCutGenerator *cg, // cut generator that has called us
127                 OsiCuts &cs,                   // cut set to be enriched
128                 int wi,
129                 int xi,
130                 CouNumber x0,           // starting point
131                 CouNumber x1,           // other bound
132                 CouNumber displacement, 
133                 bool &skip_up, 
134                 bool &skip_dn) {
135
136  CouNumber tpt,
137    rx0  = modulo (x0 + displacement, 2*M_PI),
138    rx1  = rx0 + x1 - x0,
139    base = x0 + displacement - rx0,
140    sinrx0 = sin (rx0), zero;
141
142  int
143    up   = (modulo (rx0, 2*M_PI) < M_PI) ? +1 : -1,
144    left = (x0 < x1)                     ? +1 : -1;
145
146  zero = (up>0) ? 0. : M_PI;
147
148  bool *s0, *s1;
149
150  if (up>0) {s0 = &skip_up; s1 = &skip_dn;}
151  else      {s0 = &skip_dn; s1 = &skip_up;}
152
153  if (left * (modulo (rx0, M_PI) - M_PI_2) < 0) { 
154
155    // after  flex (i.e., at \_ or /~ ) for left  bound,
156    // before flex (i.e., at _/ or ~\ ) for right bound
157
158    // out of the "belly": tangent. If on upper bay we consider the
159    // lower half-plane, and viceversa --> use -up
160    cg -> addTangent (cs, wi, xi, x0, sin (rx0), cos (rx0), -up);
161
162    // leftmost extreme to search for tangent point
163    CouNumber extr0 = .75 * M_PI * (left+1) - M_PI_2 * up; 
164
165    // in:
166    if ((left * (rx1 - M_PI * ((left - up) / 2 + 1)) <= 0) ||   // if rx1 in same "belly", or
167        (left * (rx1 - (tpt = trigNewton
168                        (rx0, extr0, extr0 + M_PI_2))) <= 0)) { // before closest leaning point
169      if (!*s1) // -> chord, if not already added in previous call
170        *s1 = (cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1,         sin (rx1), up) > 0);
171    } else     cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), up);
172  }
173  else {
174
175    // after  stationary point (i.e., _/ or ~\ ) for left bound,
176    // before stationary point (i.e., /~ or \_ ) for right bound
177 
178    //    if (left * (rx1 - left * (zero + 5*M_PI_2)) < 0) {
179    if (left * (rx1 - (4*left - up + 2) * M_PI_2) < 0) {
180      CouNumber cosrx0 = cos (rx0);
181      if (up * (sin (rx1) - sinrx0 - cosrx0 * (rx1-rx0)) < 0) 
182        // (b,sinb) below tangent --> tangent
183        cg -> addTangent (cs, wi, xi, x0, sinrx0, cosrx0, -up);
184      else {    // up: either chord or leaning plane
185        CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
186        tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
187        if (left * (rx1 - tpt) < 0) {
188          if (!*s0)
189            *s0 = cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1, sin (rx1), -up) > 0;
190        }
191        else cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), -up);
192      }
193    } else {
194      CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
195      tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
196      cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), -up);
197    }
198
199    // down: other chord or leaning plane
200    if ((left * (rx1 - (zero + M_PI)) < 0) || 
201        (left * (rx1 - (tpt = trigNewton (rx0, 
202                                          (2 +   left - up) * M_PI_2, 
203                                          (2 + 2*left - up) * M_PI_2))) < 0)) {
204      if (!*s1) 
205        *s1 = (cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1,         sin (rx1), up) > 0);
206    } else     cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), up);
207  }
208
209  return 0;
210}
Note: See TracBrowser for help on using the repository browser.