source: trunk/Couenne/src/problem/CouenneProblem.hpp @ 103

Last change on this file since 103 was 103, checked in by pbelotti, 11 years ago

fixed bug in zero-multiplicity original variables: they are set to zero in the .sol file, but they shouln't - thanks to Francois Dionne for pointing it out. split CouenneProblem?.cpp to separate constructors and reformulation from the rest. Added two new files for expression containers.

File size: 18.6 KB
Line 
1/*
2 * Name:    CouenneProblem.hpp
3 * Author:  Pietro Belotti
4 * Purpose: define the class CouenneProblem
5 *
6 * (C) Carnegie-Mellon University, 2006-09.
7 * This file is licensed under the Common Public License (CPL)
8 */
9
10#ifndef COUENNE_PROBLEM_HPP
11#define COUENNE_PROBLEM_HPP
12
13#include <vector>
14
15#include "OsiRowCut.hpp"
16#include "BonBabInfos.hpp"
17
18#include "CouenneTypes.hpp"
19#include "expression.hpp"
20#include "exprAux.hpp"
21#include "CouenneProblemElem.hpp"
22#include "CouenneObject.hpp"
23#include "CouenneJournalist.hpp"
24#include "domain.hpp"
25
26struct ASL;
27struct expr;
28
29class DepGraph;
30class CouenneCutGenerator;
31class CouenneSolverInterface;
32class quadElem;
33class LinMap;
34class QuadMap;
35
36// default tolerance for checking feasibility (and integrality) of NLP solutions
37const CouNumber feas_tolerance_default = 1e-5;
38
39
40/** Class for MINLP problems with symbolic information
41 *
42 *  It is read from an AMPL .nl file and contains variables, AMPL's
43 *  "defined variables" (aka common expressions), objective(s), and
44 *  constraints in the form of expression's. Changes throughout the
45 *  program occur in standardization.
46 */
47
48class CouenneProblem {
49
50  /// Class for storing a global cutoff for a CouenneProblem and all
51  /// its clones
52  class GlobalCutOff {
53  private:
54    GlobalCutOff(const GlobalCutOff&);
55    double cutoff_;
56  public:
57    GlobalCutOff ()         : cutoff_ (COIN_DBL_MAX) {}
58    GlobalCutOff (double c) : cutoff_ (c) {}
59    ~GlobalCutOff() {}
60    inline void setCutOff (double cutoff) {cutoff_ = cutoff;}
61    inline double getCutOff() const {return cutoff_;}
62  };
63
64  /// structure to record fixed, non-fixed, and continuous variables
65  enum fixType {UNFIXED, FIXED, CONTINUOUS};
66
67 protected:
68
69  /// problem name
70  std::string problemName_;
71
72  std::vector <exprVar           *> variables_;   ///< Variables (original, auxiliary, and defined)
73  std::vector <CouenneObjective  *> objectives_;  ///< Objectives
74  std::vector <CouenneConstraint *> constraints_; ///< Constraints
75
76  /// AMPL's common expressions (read from AMPL through structures cexps and cexps1)
77  std::vector <expression *> commonexprs_; 
78
79  mutable Domain domain_; ///< current point and bounds;
80
81  /// Expression map for comparison in standardization and to count
82  /// occurrences of an auxiliary
83  std::set <exprAux *, compExpr> *auxSet_;
84
85  /// Number of elements in the x_, lb_, ub_ arrays
86  mutable int curnvars_;
87
88  /// Number of discrete variables
89  int nIntVars_;
90
91  /// Best solution known to be loaded from file -- for testing purposes
92  mutable CouNumber *optimum_;
93
94  /// Best known objective function
95  CouNumber bestObj_;
96
97  /// Indices of variables appearing in products (used for SDP cuts)
98  int *quadIndex_;
99
100  /// Variables that have commuted to auxiliary
101  bool *commuted_;
102
103  /// numbering of variables. No variable xi with associated pi(i)
104  /// greater than pi(j) should be evaluated before variable xj
105  int *numbering_;
106
107  /// Number of "defined variables" (aka "common expressions")
108  int ndefined_;
109
110  /// Dependence (acyclic) graph: shows dependence of all auxiliary
111  /// variables on one another and on original variables. Used to
112  /// create a numbering of all variables for evaluation and bound
113  /// tightening (reverse order for implied bounds)
114  DepGraph *graph_;
115
116  /// Number of original variables
117  int nOrigVars_;
118
119  /// Number of original constraints (disregarding those that turned
120  /// into auxiliary variable definition)
121  int nOrigCons_;
122
123  /// Number of original integer variables
124  int nOrigIntVars_;
125
126  /// Pointer to a global cutoff object
127  mutable GlobalCutOff* pcutoff_;
128
129  /// flag indicating if this class is creator of global cutoff object
130  mutable bool created_pcutoff_;
131
132  bool doFBBT_;  ///< do Feasibility-based bound tightening
133  bool doOBBT_;  ///< do Optimality-based  bound tightening
134  bool doABT_;   ///< do Aggressive        bound tightening
135
136  int logObbtLev_;   ///< frequency of Optimality-based bound tightening
137  int logAbtLev_;    ///< frequency of Aggressive       bound tightening
138
139  /// SmartPointer to the Journalist
140  JnlstPtr jnlst_;
141
142  /// window around known optimum (for testing purposes)
143  CouNumber opt_window_;
144
145  /// Use quadratic expressions?
146  bool useQuadratic_;
147
148  /// feasibility tolerance (to be used in checkNLP)
149  CouNumber feas_tolerance_;
150
151  /// inverse dependence structure: for each variable x give set of
152  /// auxiliary variables (or better, their indices) whose expression
153  /// depends on x
154  std::vector <std::set <int> > dependence_;
155
156  /// vector of pointer to CouenneObjects. Used by CouenneVarObjects
157  /// when finding all objects related to (having as argument) a
158  /// single variable
159  std::vector <CouenneObject *> objects_;
160
161  /// each element is true if variable is integer and, if auxiliary,
162  /// depends on no integer
163  mutable int *integerRank_;
164
165  /// numberInRank_ [i] is the number of integer variables in rank i
166  mutable std::vector <int> numberInRank_;
167
168  /// maximum cpu time
169  double maxCpuTime_;
170
171  /// options
172  Bonmin::BabSetupBase *bonBase_;
173
174#ifdef COIN_HAS_ASL
175  /// AMPL structure pointer (temporary --- looking forward to embedding into OS...)
176  ASL *asl_;
177#endif
178
179  /// some originals may be unused due to their zero multiplicity
180  /// (that happens when they are duplicates). This array keeps track
181  /// of their indices and is sorted by evaluation order
182  int *unusedOriginalsIndices_;
183
184  /// number of unused originals
185  int nUnusedOriginals_;
186
187 public:
188
189  CouenneProblem  (ASL * = NULL,
190                   Bonmin::BabSetupBase *base = NULL,
191                   JnlstPtr jnlst = NULL);  ///< Constructor
192  CouenneProblem  (const CouenneProblem &); ///< Copy constructor
193  ~CouenneProblem ();                       ///< Destructor
194
195  /// Clone method (for use within CouenneCutGenerator::clone)
196  CouenneProblem *clone () const
197  {return new CouenneProblem (*this);}
198
199  int nObjs     () const {return (int) objectives_.   size ();} ///< Get number of objectives
200  int nCons     () const {return (int) constraints_.  size ();} ///< Get number of constraints
201  int nOrigCons () const {return nOrigCons_;}                   ///< Get number of original constraints
202
203  inline int nOrigVars    () const {return nOrigVars_;}                ///< Number of orig. variables
204  inline int nOrigIntVars () const {return nOrigIntVars_;}             ///< Number of original integers
205  inline int nIntVars     () const {return nIntVars_;}                 ///< Number of integer variables
206  inline int nVars        () const {return (int) variables_. size ();} ///< Total number of variables
207
208  /// get evaluation order index
209  inline int evalOrder (int i) const
210  {return numbering_ [i];}
211
212  /// get evaluation order vector (numbering_)
213  inline int *evalVector ()
214  {return numbering_;}
215
216  // get elements from vectors
217  inline CouenneConstraint *Con (int i) const {return constraints_ [i];} ///< i-th constraint
218  inline CouenneObjective  *Obj (int i) const {return objectives_  [i];} ///< i-th objective
219
220  /// Return pointer to i-th variable
221  inline exprVar *Var   (int i) const 
222  {return variables_ [i];}
223
224  /// Return vector of variables (symbolic representation)
225  inline std::vector <exprVar *> &Variables () 
226  {return variables_;}
227
228  /// Return pointer to set for comparisons
229  inline std::set <exprAux *, compExpr> *& AuxSet () 
230  {return auxSet_;}
231
232  /// Return pointer to dependence graph
233  inline DepGraph *getDepGraph () 
234  {return graph_;}
235
236  /// return current point & bounds
237  inline Domain *domain () const
238  {return &domain_;}
239
240  // Get and set current variable and bounds
241  inline CouNumber   &X     (int i) const {return domain_.x   (i);} ///< \f$x_i\f$
242  inline CouNumber   &Lb    (int i) const {return domain_.lb  (i);} ///< lower bound on \f$x_i\f$
243  inline CouNumber   &Ub    (int i) const {return domain_.ub  (i);} ///< upper bound on\f$x_i\f$
244
245  // get and set current variable and bounds
246  inline CouNumber  *X     () const {return domain_.();} ///< Return vector of variables
247  inline CouNumber  *Lb    () const {return domain_.lb ();} ///< Return vector of lower bounds
248  inline CouNumber  *Ub    () const {return domain_.ub ();} ///< Return vector of upper bounds
249
250  // get optimal solution and objective value
251  CouNumber  *&bestSol () const {return optimum_;} ///< Best known solution (read from file)
252  CouNumber    bestObj () const {return bestObj_;} ///< Objective of best known solution
253
254  /// Get vector of commuted variables
255  bool *&Commuted () 
256  {return commuted_;}
257
258  /// Add (non linear) objective function
259  void addObjective     (expression *, const std::string &);
260
261  // Add (non linear) "=", ">=", "<=", and range constraints
262  void addEQConstraint  (expression *, expression *); ///< Add equality constraint \f$ h(x) = b\f$
263  void addGEConstraint  (expression *, expression *); ///< Add \f$\ge\f$ constraint, \f$h(x)\ge b\f$
264  void addLEConstraint  (expression *, expression *); ///< Add \f$\le\f$ constraint, \f$h(x)\le b\f$
265  void addRNGConstraint (expression *, expression *, 
266                         expression *);               ///< Add range constraint, \f$a\le h(x)\le b\f$
267
268  /// Add original variable.
269  ///
270  /// @param isint if true, this variable is integer, otherwise it is
271  /// continuous
272  expression *addVariable (bool isint = false, Domain *d = NULL);
273
274  /// Add auxiliary variable and associate it with expression given as
275  /// argument (used in standardization)
276  exprAux *addAuxiliary (expression *);
277
278  /// preprocess problem in order to extract linear relaxations etc.
279  void reformulate ();
280
281  /// Break problem's nonlinear constraints in simple expressions to
282  /// be convexified later. Return true if problem looks feasible,
283  /// false if proven infeasible.
284  bool standardize ();
285
286  /// Display current representation of problem: objective, linear and
287  /// nonlinear constraints, and auxiliary variables.
288  void print (std::ostream & = std::cout);
289
290#ifdef COIN_HAS_ASL
291  /// Read problem from .nl file using the Ampl Solver Library (ASL)
292  int readnl (const struct ASL *);
293
294  /// Generate a Couenne expression from an ASL expression
295  expression *nl2e (struct expr *, const ASL *asl);
296#endif
297
298  // bound tightening parameters
299  bool doFBBT () const {return doFBBT_;} ///< shall we do Feasibility Based Bound Tightening?
300  bool doOBBT () const {return doOBBT_;} ///< shall we do Optimality  Based Bound Tightening?
301  bool doABT  () const {return doABT_;}  ///< shall we do Aggressive        Bound Tightening?
302
303  int  logObbtLev () const {return logObbtLev_;} ///< How often shall we do OBBT?
304  int  logAbtLev  () const {return logAbtLev_;}  ///< How often shall we do ABT?
305
306  /// Write nonlinear problem to a .mod file (with lots of defined
307  /// variables)
308  ///
309  /// @param fname Name of the .mod file to be written
310  ///
311  /// @param aux controls the use of auxiliaries. If true, a problem
312  /// is written with auxiliary variables written with their
313  /// associated expression, i.e. \f$w_i = h_i(x,y,w)\f$ and bounds
314  /// \f$l_i \le w_i \le u_i\f$, while if false these constraints are
315  /// written in the form \f$l_i \le h_i (x,y) \le u_i\f$.
316  ///
317  /// Note: if used before standardization, writes original AMPL formulation
318  void writeAMPL (const std::string &fname, bool aux);
319
320  /// Write nonlinear problem to a .gms file
321  ///
322  /// @param fname Name of the .gams file to be written.
323  void writeGAMS (const std::string &fname);
324
325  /// Initialize auxiliary variables and their bounds from original
326  /// variables
327  //void initAuxs (const CouNumber *, const CouNumber *, const CouNumber *);
328  void initAuxs () const;
329
330  /// Get auxiliary variables from original variables
331  void getAuxs (CouNumber *) const;
332
333  /// tighten bounds using propagation, implied bounds and reduced costs
334  bool boundTightening (t_chg_bounds *, 
335                        Bonmin::BabInfo * = NULL) const;
336
337  /// core of the bound tightening procedure
338  bool btCore (t_chg_bounds *chg_bds) const;
339
340  /// Optimality Based Bound Tightening
341  int obbt (const CouenneCutGenerator *cg,
342            const OsiSolverInterface &csi,
343            OsiCuts &cs,
344            const CglTreeInfo &info,
345            Bonmin::BabInfo * babInfo,
346            t_chg_bounds *chg_bds);
347
348  /// aggressive bound tightening. Fake bounds in order to cut
349  /// portions of the solution space by fathoming on
350  /// bounds/infeasibility
351  bool aggressiveBT (Bonmin::OsiTMINLPInterface *nlp,
352                     t_chg_bounds *, 
353                     Bonmin::BabInfo * = NULL) const;
354
355  /// procedure to strengthen variable bounds. Return false if problem
356  /// turns out to be infeasible with given bounds, true otherwise.
357  int redCostBT (const OsiSolverInterface *psi,
358                 t_chg_bounds *chg_bds) const;
359
360  /// "Forward" bound tightening, that is, propagate bound of variable
361  /// \f$x\f$ in an expression \f$w = f(x)\f$ to the bounds of \f$w\f$.
362  int tightenBounds (t_chg_bounds *) const;
363
364  /// "Backward" bound tightening, aka implied bounds.
365  int impliedBounds (t_chg_bounds *) const;
366
367  /// Look for quadratic terms to be used with SDP cuts
368  void fillQuadIndices ();
369
370  /// Fill vector with coefficients of objective function
371  void fillObjCoeff (double *&);
372
373  /// Replace all occurrences of original variable with new aux given
374  /// as argument
375  void auxiliarize (exprVar *, exprVar * = NULL);
376
377  /// Set cutoff
378  void setCutOff (CouNumber cutoff) const;
379
380  /// Set cutoff
381  CouNumber getCutOff () const
382  {return pcutoff_ -> getCutOff ();}
383
384  /// Make cutoff known to the problem
385  void installCutOff () const;
386
387  /// Provide Journalist
388  ConstJnlstPtr Jnlst() const 
389  {return ConstPtr(jnlst_);}
390
391  /// Check if solution is MINLP feasible
392  bool checkNLP (const double *solution, double &obj, bool recompute = false) const;
393
394  /// generate integer NLP point Y starting from fractional solution
395  /// using bound tightening
396  int getIntegerCandidate (const double *xFrac, double *xInt, double *lb, double *ub) const;
397
398  /// Read best known solution from file given in argument
399  bool readOptimum (std::string *fname = NULL);
400
401  /// Read cutoff value (for testing purposes)
402  void readCutoff (const std::string &fname);
403
404  /// Add list of options to be read from file
405  static void registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions);
406
407  /// standardization of linear exprOp's
408  exprAux *linStandardize (bool addAux, 
409                           CouNumber c0, 
410                           LinMap  &lmap,
411                           QuadMap &qmap);
412
413  /// split a constraint w - f(x) = c into w's index (it is returned)
414  /// and rest = f(x) + c
415  int splitAux (CouNumber, expression *, expression *&, bool *);
416
417  /// translates pair (indices, coefficients) into vector with pointers to variables
418  void indcoe2vector (int *indexL,
419                      CouNumber *coeff,
420                      std::vector <std::pair <exprVar *, CouNumber> > &lcoeff);
421
422  /// translates triplet (indicesI, indicesJ, coefficients) into vector with pointers to variables
423  void indcoe2vector (int *indexI,
424                      int *indexJ,
425                      CouNumber *coeff,
426                      std::vector <quadElem> &qcoeff);
427
428  /// given (expression *) element of sum, returns (coe,ind0,ind1)
429  /// depending on element:
430  ///
431  /// 1) a * x_i ^ 2   ---> (a,i,?)   return COU_EXPRPOW
432  /// 2) a * x_i       ---> (a,i,?)   return COU_EXPRVAR
433  /// 3) a * x_i * x_j ---> (a,i,j)   return COU_EXPRMUL
434  /// 4) a             ---> (a,?,?)   return COU_EXPRCONST
435  ///
436  /// x_i and/or x_j may come from standardizing other (linear or
437  /// quadratic operator) sub-expressions
438  void decomposeTerm (expression *term,
439                      CouNumber initCoe,
440                      CouNumber &c0,
441                      LinMap  &lmap,
442                      QuadMap &qmap);
443
444  /// return problem name
445  const std::string &problemName () const
446  {return problemName_;}
447
448  /// return inverse dependence structure
449  const std::vector <std::set <int> > &Dependence () const
450  {return dependence_;}
451
452  /// return object vector
453  const std::vector <CouenneObject *> &Objects () const
454  {return objects_;}
455
456  /// find SOS constraints in problem
457  int findSOS (OsiSolverInterface *solver, OsiObject ** objects);
458
459  /// set maximum CPU time
460  inline void setMaxCpuTime (double time)
461  {maxCpuTime_ = time;}
462
463  /// return maximum CPU time
464  inline double getMaxCpuTime () const
465  {return maxCpuTime_;}
466
467  /// save CouenneBase
468  inline void setBase (Bonmin::BabSetupBase *base) {
469    bonBase_ = base;
470    jnlst_   = base -> journalist ();
471  }
472
473  /// Some originals may be unused due to their zero multiplicity
474  /// (that happens when they are duplicates). This procedure creates
475  /// a structure for quickly checking and restoring their value after
476  /// solving.
477  void createUnusedOriginals ();
478
479  /// Some originals may be unused due to their zero multiplicity (that
480  /// happens when they are duplicates). This procedure restores their
481  /// value after solving
482  void restoreUnusedOriginals (CouNumber * = NULL) const;
483
484  /// return indices of neglected redundant variables
485  int *unusedOriginalsIndices () {return unusedOriginalsIndices_;}
486
487  /// number of unused originals
488  int nUnusedOriginals ()        {return nUnusedOriginals_;}
489
490protected:
491
492  /// single fake tightening. Return
493  ///
494  /// -1   if infeasible
495  ///  0   if no improvement
496  /// +1   if improved
497  int fake_tighten (char direction,  ///< 0: left, 1: right
498                    int index,       ///< index of the variable tested
499                    const double *X, ///< point round which tightening is done
500                    CouNumber *olb,  ///< cur. lower bound
501                    CouNumber *oub,  ///< cur. upper bound
502                    t_chg_bounds *chg_bds,
503                    t_chg_bounds *f_chg) const;
504
505  /// Optimality Based Bound Tightening -- inner loop
506  int obbtInner (CouenneSolverInterface *, 
507                 OsiCuts &,
508                 t_chg_bounds *, 
509                 Bonmin::BabInfo *) const;
510
511  int obbt_iter (CouenneSolverInterface *csi, 
512                 t_chg_bounds *chg_bds, 
513                 const CoinWarmStart *warmstart, 
514                 Bonmin::BabInfo *babInfo,
515                 double *objcoe,
516                 int sense, 
517                 int index) const;
518
519  int call_iter (CouenneSolverInterface *csi, 
520                 t_chg_bounds *chg_bds, 
521                 const CoinWarmStart *warmstart, 
522                 Bonmin::BabInfo *babInfo,
523                 double *objcoe,
524                 enum nodeType type,
525                 int sense) const;
526
527  /// analyze sparsity of potential exprQuad/exprGroup and change
528  /// linear/quadratic maps accordingly, if necessary by adding new
529  /// auxiliary variables and including them in the linear map
530  void analyzeSparsity (CouNumber, 
531                        LinMap &,
532                        QuadMap &);
533
534  /// re-organizes multiplication and stores indices (and exponents) of
535  /// its variables
536  void flattenMul (expression *mul, 
537                   CouNumber &coe, 
538                   std::map <int, CouNumber> &indices);
539
540  /// clear all spurious variables pointers not referring to the variables_ vector
541  void realign ();
542
543  /// fill dependence_ structure
544  void fillDependence (Bonmin::BabSetupBase *base);
545
546  /// fill freeIntegers_ array
547  void fillIntegerRank () const;
548
549  /// Test fixing of an integer variable (used in getIntegerCandidate())
550  int testIntFix (int index, 
551                  CouNumber xFrac, 
552                  enum fixType *fixed,
553                  CouNumber *xInt,
554                  CouNumber *dualL, CouNumber *dualR,
555                  CouNumber *olb,   CouNumber *oub,
556                  bool patient) const;
557};
558
559#endif
Note: See TracBrowser for help on using the repository browser.