simplex.c
Go to the documentation of this file.
1 #include "simplex.h"
2 
3 #include "defines.h"
4 #include "random.h"
5 
6 
19 void SimplexExpandBounds(unsigned int eq_type,Tinterval *b)
20 {
21  double v1,v2;
22 
23  v1=LowerLimit(b);
24  v2=UpperLimit(b);
25 
26  if (eq_type==LEQ)
27  v1=-INF;
28  else
29  {
30  if (eq_type==GEQ)
31  v2=INF;
32  }
33 
34  NewInterval(v1,v2,b);
35 }
36 
37 /*
38  Sets new ranges for the variables in the simplex
39 */
41 {
42  unsigned int i;
43  Tinterval *is;
44  unsigned int nv;
45 
46  #if (_DEBUG>4)
47  printf(" Setting bound for columns (variable):\n");
48  #endif
49 
50  nv=GetBoxNIntervals(b);
51  is=GetBoxIntervals(b);
52  for(i=0;i<nv;i++)
53  {
54  SimplexSetColBounds(i,&(is[i]),lp);
55 
56  #if (_DEBUG>4)
57  printf(" %u: ",i);
58  PrintInterval(stdout,&(is[i]));
59  if (IntervalSize(&(is[i]))<INF)
60  printf(" -> %g\n",IntervalSize(&(is[i])));
61  else
62  printf(" -> +inf\n");
63  #endif
64  }
65 }
66 
67 boolean SimplexAddNewConstraint(double epsilon,
68  unsigned int safeSimplex,
70  unsigned int eq_type,
71  Tinterval *is,
72  TSimplex *s)
73 {
74  boolean full;
75 
76  full=TRUE;
77 
79  {
80  TLinearConstraint lcInt;
81  unsigned int n;
82  Tinterval bounds;
83  double es,ic;
84 
85  /*We avoid any modification in the given linear constraint. Modification
86  can produce undesired side-effects out of this method*/
87  CopyLinearConstraint(&lcInt,lc);
88 
89  #if (_DEBUG>4)
90  printf(" Adding new constraint:\n ");
91  PrintLinearConstraint(stdout,TRUE,NULL,&lcInt);
92  #endif
93 
94  /* Small ranges for the variables can produce simplex inestabilities.
95  To improve numerical stability we can subsume those small ranges in
96  the right-hand side of the equation.
97  For inequalities this operation is done for safeSimplex>=6.
98  For equalities this is only done for safeSimplex>=5. This is so to
99  minimize the risk of converting an equality into an inequality.
100 
101  For GLPK it is compulsatory to use this. For the other simplex engines
102  a smaller threshold could be used. */
103 
104  if ((safeSimplex>=6)||
105  ((safeSimplex>=5)&&(GetLinearConstraintErrorSize(&lcInt)>0)))
106  CleanLinearConstraint((_SIMPLEX_ENGINE==_GLPK?epsilon:ZERO),is,&lcInt);
107 
108  #if (_DEBUG>4)
109  printf(" Constraint after cleaning:\n ");
110  PrintLinearConstraint(stdout,TRUE,NULL,&lcInt);
111  #endif
112 
114  /* If the input linear constraint involves at least one variable...*/
115  if (n>0)
116  {
117  /* linear constraint with very narrow range are converted to inequalities
118  (randomly selecting LEQ or GEQ). This is a conservative way to avoid
119  numerical issues. */
120  GetLinearConstraintError(&bounds,&lcInt);
121  SimplexExpandBounds(eq_type,&bounds);
122  SetLinearConstraintError(&bounds,&lcInt);
123 
124  #if (_DEBUG>4)
125  printf(" Constraint after expanding bounds:\n ");
126  PrintLinearConstraint(stdout,TRUE,NULL,&lcInt);
127  #endif
128 
129  if (!SimplifyLinearConstraint(&full,is,&lcInt))
130  {
131  /*******************************************************************/
132  /* Extremely small ranges for equalities are not correctly processed
133  by the simplex engines (they convert the tiny ranges to a point and
134  then solutions might be lost)
135 
136  Those small ranges are produced in problems where we use a tiny
137  epsilon or due to the CleanLinearConstraint process (possibly) used
138  above.
139  */
140  if (eq_type==EQU)
141  {
142  es=IntervalSize(&bounds);
143  if(safeSimplex>=4)
144  {
145  if (es>0)
146  {
147  if (es<=ZERO)
148  {
149  /*This is an equality that somehow has been converted in
150  a double bounded constraint with very tiny range (or error).
151  We convert it back to an equality.
152  Note that with this we can lose solutions.
153  Use safeSimplex<4 to avoid this process.*/
154  ic=IntervalCenter(&bounds);
155  NewInterval(ic,ic,&bounds);
156  SetLinearConstraintError(&bounds,&lcInt);
157  }
158  else
159  {
160  if ((safeSimplex>=7)&&(es<epsilon))
161  //if ((safeSimplex>=7)&&(es<(_SIMPLEX_ENGINE==_GLPK?epsilon:ZERO)))
162  {
163  /*For double bounded constraints with very tiny
164  error, we use only one of the constraints >= or <=
165  selected randomly. This is numerically safer.
166  However, if the user selected an epsilon below 1e-6
167  we assume she wants to be in the wild side and
168  that constraint errors can be arbitrarily small.
169  */
170  SimplexExpandBounds((randomDouble()<0.5?LEQ:GEQ),&bounds);
171  SetLinearConstraintError(&bounds,&lcInt);
172  }
173  }
174  }
175  }
176  /*
177  else
178  {
179  if (es<epsilon)
180  {
181  SimplexExpandBounds((randomDouble()<0.5?LEQ:GEQ),&bounds);
182  SetLinearConstraintError(&bounds,&lcInt);
183  //Tinterval kk;
184 
185  //NewInterval(-epsilon/2.0,epsilon/2.0,&kk);
186  //IntervalAdd(&bounds,&kk,&bounds);
187  //SetLinearConstraintError(&bounds,&lcInt);
188  }
189  }
190  */
191  }
192  /*******************************************************************/
193 
194  SimplexAddNewConstraintRaw(&lcInt,s);
195  }
196  #if (_DEBUG>4)
197  else
198  printf(" Empty simplified equation\n");
199  #endif
200  }
201  #if (_DEBUG>4)
202  else
203  printf(" Empty input equation\n");
204  #endif
205 
206  DeleteLinearConstraint(&lcInt);
207  }
208  #if (_DEBUG>4)
209  else
210  printf(" Unbounded linear constraint. No constraint.\n");
211  #endif
212  return(full);
213 }
214 
215 double SimplexGetOptimalValue(unsigned int safeSimplex,double m,Tbox *x,TSimplex *s)
216 {
217  double o_safe;
218 
219  if (safeSimplex==0)
220  o_safe=SimplexGetOptimalValueRaw(s);
221  else
222  {
223  /*
224  LP safe bounds taken from
225  "Safe bounds in linear and mixed-integer linear programming"
226  Arnold Neumaier and Oleg Shcherbina, 2003
227  */
228 
229  double *lambda,*obj_c;
230  double o,rl,ru;
231  unsigned int n,m;
232  unsigned int i,j,k;;
233  Tinterval r,obj,a;
235  Tinterval bounds;
236  double ct_l,ct_u;
237  TLinearConstraint s_obj;
238 
239  m=SimplexNRows(s);
240  n=SimplexNColumns(s);
241 
242  NEW(lambda,m,double);
243  for(j=0;j<m;j++)
244  lambda[j]=SimplexGetRowDual(j,s);
245 
246  NEW(obj_c,n,double);
247  for(i=0;i<n;i++)
248  obj_c[i]=0.0;
251  for(i=0;i<k;i++)
253  DeleteLinearConstraint(&s_obj);
254 
255  /* o_safe = inf{ (\sum_j^m lambda[j]*b[j]) - (sum_i^n r[i]*x[i])}
256  with
257  r[i] = - c[i] + sum_j^m A(j,i)*lambda[j]
258  Note that only the lower limit is necessary, thus
259  o_safe = inf{\sum_j^m lambda[j]*b[j]} -sup{sum_i^n r[i]*x[i]}
260  */
261 
262  o_safe=0.0; /*= inf{\sum_j^m lambda[j]*b[j]} */
263  ROUNDDOWN;
264  for(j=0;j<m;j++)
265  {
266  SimplexGetRowBounds(j,&bounds,s);
267 
268  ct_l=LowerLimit(&bounds);
269  ct_u=UpperLimit(&bounds);
270 
271  if ((ct_u==INF)||(lambda[j]>0))
272  o_safe+=(lambda[j]*ct_l);
273  else
274  o_safe+=(lambda[j]*ct_u);
275  }
276  ROUNDNEAR;
277 
278  NewInterval(0,0,&obj); /*obj=sum_i^n r[i]*x[i]*/
279 
280  for(i=0;i<n;i++)
281  {
282  /*Compute r[i] = -c[i] + sum_j^m A(j,i)*lambda[j] */
283  SimplexGetColConstraint(i,&lc,s);
285 
286  ROUNDDOWN;
287  rl=-obj_c[i];
288  for(j=0;j<k;j++)
290 
291  ROUNDUP;
292  ru=-obj_c[i];
293  for(j=0;j<k;j++)
295 
296  ROUNDNEAR;
297  NewInterval(rl,ru,&r);
298 
299  /*Compute r[i]*x[i]*/
300  IntervalProduct(&r,GetBoxInterval(i,x),&a);
301 
302  /*Add to obj*/
303  IntervalAdd(&obj,&a,&obj);
304 
306  }
307 
308  free(lambda);
309  free(obj_c);
310 
311  ROUNDDOWN;
312  o_safe=o_safe-UpperLimit(&obj);
313  ROUNDNEAR;
314 
316 
317  o_safe=(o_safe<o?o_safe:o);
318  }
319 
320  return((o_safe<m?m:o_safe));
321 }
322 
323 
324 /*
325  Reduces a box along a given range. The reduction is based on minimizing
326  the simplex problem with a objective function that only affects the given
327  varible.
328 */
329 unsigned int ReduceRange(double epsilon,unsigned int safeSimplex,
330  unsigned int nv,Tbox *b,TSimplex *lp)
331 {
332  unsigned int j;
333  unsigned int r; /*Simplex result*/
334  boolean empty;
335  boolean err; /*Simplex error*/
336  Tinterval new_range;
337  Tinterval *i_in;
338  double s_in;
339  double o;
340  TLinearConstraint obj;
341 
342  empty=FALSE;
343  err=FALSE;
344 
345  /* Add the most up to date bounds for the system/dummy variables to the simplex */
346  SetSimplexBounds(b,lp);
347 
348  #if (_SIMPLEX_ENGINE!=_CLP) /* in CLP Reset Simplex has no effect */
349  if (safeSimplex>=2)
350  {
351  /* This improves the numerical stability but slows down the process */
352  ResetSimplex(lp);
353  }
354  #endif
355 
356  i_in=GetBoxInterval(nv,b); /* Range to be reduced */
357 
358  s_in=IntervalSize(i_in);
359 
360  #if (_DEBUG>4)
361  printf(" Reducing interval %u via simplex\n",nv);
362  printf(" Original Interval ");
363  PrintInterval(stdout,i_in);
364  if (s_in>=INF) printf(" s: inf\n");
365  else printf(" s:%g\n",s_in);
366  #endif
367 
368  CopyInterval(&new_range,i_in); /* Default initialization used in case
369  of ERROR_IN_PROCESS*/
370  if (s_in>=epsilon)
371  {
372  InitLinearConstraint(&obj);
373 
374  for(j=0;((!empty)&&(!err)&&(j<2));j++) /*j=0 minimize j=1 maximize*/
375  {
376  ResetLinearConstraint(&obj);
377  if (j==0)
378  {
379  AddTerm2LinearConstraint(nv,+1.0,&obj);
380 
381  #if (_DEBUG>4)
382  printf("Determining lower bound for variable %u -> \n",nv);
383  #endif
384  }
385  else
386  {
387  AddTerm2LinearConstraint(nv,-1.0,&obj);
388 
389  #if (_DEBUG>4)
390  printf("Determining upper bound for variable %u -> \n",nv);
391  #endif
392  }
393 
395 
396  #if (_SIMPLEX_ENGINE!=_CLP) /* in CLP Reset Simplex has no effect */
397  if ((safeSimplex>=3)&&(j==1))
398  {
399  /* This improves (even more) the numerical stability but slows down (even more)
400  the process (See the use of ResetSimplex in ReduceBox) */
401  ResetSimplex(lp);
402  }
403  #endif
404 
405  /*Solve the simplex problem*/
406  #if (_DEBUG>4)
407  printf("==================================================================================\n");
408  #endif
409  r=SimplexSolve(lp);
410 
411  #if (_SIMPLEX_ENGINE!=_CLP) /* in CLP Reset Simplex has no effect */
412  if ((r==EMPTY_BOX)||
413  (r==ERROR_IN_PROCESS)||
414  ((j==1)&&(r==REDUCED_BOX)&&(-o<LowerLimit(&new_range))))
415  {
416  /* A second chance but with clean internal
417  structures. This minimizes que change of wrongly
418  classifying a box as empty */
419  ResetSimplex(lp);
420  r=SimplexSolve(lp);
421  }
422  #endif
423 
424  if (r==REDUCED_BOX)
425  {
426  if (IntervalSize(i_in)<INF)
427  o=SimplexGetOptimalValue(safeSimplex,(j==0?LowerLimit(i_in):-UpperLimit(i_in)),b,lp);
428  else
430  }
431 
432  #if (_DEBUG>4)
433  printf("==================================================================================\n");
434  #endif
435 
436  switch(r)
437  {
438  case REDUCED_BOX:
439  if (j==0) /*min*/
440  {
441  /*a new lower: we only use the new lower limit if it is
442  higher than the current one. Lower lower limits can
443  appear due to numerical errors*/
444  UpdateLowerLimit(o,&new_range);
445 
446  #if (_DEBUG>4)
447  printf(" The new lower limit is %g\n",LowerLimit(&new_range));
448  #endif
449  }
450  else /*max*/
451  {
452  /*a new upper: we only use the new upper limit if it is
453  lower than the current one. Higer upper limits can
454  appear due to numerical errors*/
455  UpdateUpperLimit(-o,&new_range);
456 
457  #if (_DEBUG>4)
458  printf(" The new upper limit is %g\n",UpperLimit(&new_range));
459  fflush(stdout);
460  #endif
461  }
462  break;
463 
464  #if (_DEBUG>4)
465  case UNBOUNDED_BOX:
466  printf(" The range can not be bounded\n");
467  break;
468  #endif
469 
470  case EMPTY_BOX:
471  empty=TRUE;
472  break;
473 
474  case ERROR_IN_PROCESS:
475  #if (_SIMPLEX_ENGINE!=_CLP) /* in CLP Reset Simplex has no effect */
476  ResetSimplex(lp);
477  #endif
478  err=TRUE;
479  break;
480 
481  default:
482  Error("Unknow Error in ReduceRange");
483  }
484  } /*min/max*/
485 
487 
488 
489  #if (_DEBUG>4)
490  if (!empty)
491  printf(" Interval out of simplex [%g %g] s:%g\n",
492  LowerLimit(&new_range),
493  UpperLimit(&new_range),
494  IntervalSize(&new_range));
495  else
496  printf(" Empty Interval out of simplex\n");
497  #endif
498 
499 
500  if ((!empty)&&(!err))
501  {
502  /*
503  In some cases the simplex return empty intervals. This occurs due to
504  numerical inestabilities in the simplex tableau monomialization process.
505  We treat this as an error -> the box is bisected and a brand new
506  simplex-based reduction is triggered
507  */
508  if (EmptyInterval(&new_range))
509  err=TRUE;
510  else
511  CopyInterval(i_in,&new_range);
512  }
513  }
514  #if (_DEBUG>4)
515  else
516  printf(" Input interval is too small to be reduced any more\n");
517  #endif
518 
519  if (empty)
520  return(EMPTY_BOX);
521  else
522  {
523  if (err)
524  return(ERROR_IN_PROCESS);
525  else
526  return(REDUCED_BOX);
527  }
528 }
unsigned int SimplexNRows(TSimplex *s)
Gets the number of rows (i.e., constraints) of the simplex structure.
Definition: simplex_clp.c:107
void SimplexGetColConstraint(unsigned int ncol, TLinearConstraint *lc, TSimplex *s)
Gets a column from the simplex in the form of a linear constraint.
Definition: simplex_clp.c:133
void DeleteLinearConstraint(TLinearConstraint *lc)
Destructor.
Tinterval * GetBoxInterval(unsigned int n, Tbox *b)
Returns a pointer to one of the intervals defining the box.
Definition: box.c:270
boolean SimplifyLinearConstraint(boolean *full, Tinterval *is, TLinearConstraint *lc)
Apply linear constraints to reduce the ranges of the problem variables.
#define ERROR_IN_PROCESS
One of the possible results of reducing a box.
Definition: box.h:32
#define FALSE
FALSE.
Definition: boolean.h:30
void SetSimplexBounds(Tbox *b, TSimplex *lp)
Sets the columns bounds for the simplex.
Definition: simplex.c:40
#define ROUNDDOWN
Sets the floating point operations in rounding down mode.
Definition: defines.h:219
#define UNBOUNDED_BOX
One of the possible results of reducing a box.
Definition: box.h:45
#define NEW(_var, _n, _type)
Allocates memory space.
Definition: defines.h:385
void SimplexSetColBounds(unsigned int ncol, Tinterval *i, TSimplex *s)
Sets the bounds for a given column (i.e., variable).
Definition: simplex_clp.c:113
A linear constraint with an associated error.
void IntervalAdd(Tinterval *i1, Tinterval *i2, Tinterval *i_out)
Addition of two intervals.
Definition: interval.c:418
#define EQU
In a Tequation, the equation relational operator is equal.
Definition: equation.h:201
void InitLinearConstraint(TLinearConstraint *lc)
Constructor.
#define EMPTY_BOX
One of the possible results of reducing a box.
Definition: box.h:25
boolean SimplexAddNewConstraint(double epsilon, unsigned int safeSimplex, TLinearConstraint *lc, unsigned int eq_type, Tinterval *is, TSimplex *s)
Adds a row (i.e., a constraint) to the simplex.
Definition: simplex.c:67
unsigned int GetLinearConstraintVariable(unsigned int i, TLinearConstraint *lc)
Gets the a particular variable index.
void PrintLinearConstraint(FILE *f, boolean eq, char **varName, TLinearConstraint *lc)
Prints a linear constraint.
#define TRUE
TRUE.
Definition: boolean.h:21
#define LEQ
In a Tequation, the equation relational operator is less equal.
Definition: equation.h:195
void Error(const char *s)
General error function.
Definition: error.c:80
A simplex tableau structure.
Definition: simplex.h:73
#define ZERO
Floating point operations giving a value below this constant (in absolute value) are considered 0...
Definition: defines.h:37
void AddTerm2LinearConstraint(unsigned int ind, double val, TLinearConstraint *lc)
Adds a scaled variable to the linear constraint.
void CopyInterval(Tinterval *i_dst, Tinterval *i_org)
Copy constructor.
Definition: interval.c:59
double SimplexGetOptimalValue(unsigned int safeSimplex, double m, Tbox *x, TSimplex *s)
Returns the optimal value determined by the simplex corrected to compensate for possible rounding eff...
Definition: simplex.c:215
#define GEQ
In a Tequation, the equation relational operator is great equal.
Definition: equation.h:189
unsigned int GetBoxNIntervals(Tbox *b)
Box dimensionality.
Definition: box.c:992
double SimplexGetRowDual(unsigned int nrow, TSimplex *s)
Gets a row dual value after solving the simplex.
Definition: simplex_clp.c:250
double randomDouble()
Returns a random double in the [0,1] interval.
Definition: random.c:33
void ResetLinearConstraint(TLinearConstraint *lc)
Resets a linear constraint.
unsigned int SimplexSolve(TSimplex *s)
Determines an optimal value.
Definition: simplex_clp.c:332
void ResetSimplex(TSimplex *s)
Resets the simplex structure.
Definition: simplex_clp.c:96
void CopyLinearConstraint(TLinearConstraint *lc_dst, TLinearConstraint *lc_src)
Copy constructor.
void UpdateLowerLimit(double c, Tinterval *i)
Updates the lower limit.
Definition: interval.c:247
double LowerLimit(Tinterval *i)
Gets the lower limit.
Definition: interval.c:79
void SimplexExpandBounds(unsigned int eq_type, Tinterval *b)
Expands an interval according to the equation type.
Definition: simplex.c:19
Definitions of constants and macros used in several parts of the cuik library.
unsigned int SimplexNColumns(TSimplex *s)
Gets the number of columns (i.e., variables) of the simplex structure.
Definition: simplex_clp.c:102
boolean BoundedLinearConstraint(TLinearConstraint *lc)
Test if the constraint is bounded.
#define ROUNDNEAR
Sets the floating point operations in round near mode.
Definition: defines.h:225
void SimplexAddNewConstraintRaw(TLinearConstraint *lc, TSimplex *s)
Adds a row (i.e., a constraint) to the simplex.
Definition: simplex_clp.c:259
double SimplexGetOptimalValueRaw(TSimplex *s)
Gets the optimal value after optimizing the problem.
Definition: simplex_clp.c:327
double UpperLimit(Tinterval *i)
Gets the uppser limit.
Definition: interval.c:87
void SimplexGetRowBounds(unsigned int nrow, Tinterval *i, TSimplex *s)
Gets the bounds for a given row (i.e., constraint).
Definition: simplex_clp.c:193
unsigned int GetNumTermsInLinearConstraint(TLinearConstraint *lc)
Number of variables in a linear constraint.
A box.
Definition: box.h:83
void GetLinearConstraintError(Tinterval *error, TLinearConstraint *lc)
Gets the right-hand side interval for the linear constraint.
void SimplexSetOptimizationFunction(TLinearConstraint *obj, TSimplex *s)
Sets a new objective function.
Definition: simplex_clp.c:286
Definition of the TSimplex type and the associated functions.
void PrintInterval(FILE *f, Tinterval *i)
Prints an interval.
Definition: interval.c:1001
double IntervalCenter(Tinterval *i)
Gets the interval center.
Definition: interval.c:129
#define ROUNDUP
Sets the floating point operations in rounding up mode.
Definition: defines.h:213
boolean EmptyInterval(Tinterval *i)
Checks if a given interval is empty.
Definition: interval.c:335
Tinterval * GetBoxIntervals(Tbox *b)
Returns a pointer to the array of intervals defining the box.
Definition: box.c:284
void SetLinearConstraintError(Tinterval *error, TLinearConstraint *lc)
Sets a new righ-hand side error of the linear constraint.
void IntervalProduct(Tinterval *i1, Tinterval *i2, Tinterval *i_out)
Product of two intervals.
Definition: interval.c:384
void NewInterval(double lower, double upper, Tinterval *i)
Constructor.
Definition: interval.c:47
#define INF
Infinite.
Definition: defines.h:70
double GetLinearConstraintCoefficient(unsigned int i, TLinearConstraint *lc)
Gets the a particular linear constraint coefficient.
unsigned int ReduceRange(double epsilon, unsigned int safeSimplex, unsigned int nv, Tbox *b, TSimplex *lp)
Reduces a variable range using the simplex.
Definition: simplex.c:329
#define REDUCED_BOX
One of the possible results of reducing a box.
Definition: box.h:38
Definition of basic randomization functions.
void SimplexGetOptimizationFunction(TLinearConstraint *obj, TSimplex *s)
Gets a current objective function.
Definition: simplex_clp.c:312
Defines a interval.
Definition: interval.h:33
void UpdateUpperLimit(double c, Tinterval *i)
Updates the upper limit.
Definition: interval.c:253
double IntervalSize(Tinterval *i)
Gets the uppser limit.
Definition: interval.c:96
void CleanLinearConstraint(double epsilon, Tinterval *is, TLinearConstraint *lc)
Removes terms in the linear constraint that give too small ranges.
double GetLinearConstraintErrorSize(TLinearConstraint *lc)
Gets the size of the right-hand side interval for the linear constraint.