solve.c

利用c语言编写

    drow[0] = 1;
    btran(lp, drow, lp->epsel);
    for(i = 1; i <= lp->columns; i++) {
      varnr = lp->rows + i;
      if(!lp->basis[varnr]) {
	f = 0;
	for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++)
	  f += drow[lp->mat[j].row_nr] * lp->mat[j].value;
	drow[varnr] = f;
      }
    }
    for(i = 1; i <= lp->sum; i++)
      my_round(drow[i], lp->epsd);

    /* original (unscaled) objective function */
    get_row(lp, 0, OrigObj);

    infinite=lp->infinite;
    epsel=lp->epsel;
    for(i = 1; i <= lp->columns; i++) {
     from=-infinite;
     till= infinite;
     varnr = lp->rows + i;
     if (!lp->basis[varnr]) {
      /* only the coeff of the objective function of column i changes. */
      a=drow[varnr];
      if (lp->scaling_used)
       a/=(lp->scale[varnr]*lp->scale[0]);
      if (lp->upbo[varnr]==0.0) /* ignore, because this case doesn't results in further iterations */ ;
      else if (lp->lower[varnr]) from=OrigObj[i]-a; /* less than this value gives further iterations */
      else till=OrigObj[i]-a; /* bigger than this value gives further iterations */
     }
     else {
      /* all the coeff of the objective function change. Search the minimal change needed for further iterations */
      for (row_nr=1;(row_nr<=lp->rows) && (lp->bas[row_nr]!=varnr);row_nr++); /* search on which row the variable exists in the basis */
      if (row_nr<=lp->rows) {     /* safety test; should always be found ... */
       setpivrow(lp,row_nr,prow); /* construct one row of the tableau */
       min1=infinite;
       min2=infinite;
       for (l=1;l<=lp->sum;l++)   /* search for the column(s) which first results in further iterations */
	if ((!lp->basis[l]) && (lp->upbo[l]>0.0) && (my_abs(prow[l])>epsel) && (drow[l]*(lp->lower[l] ? -1 : 1)<epsel)) {
	 a=my_abs(drow[l]/prow[l]);
	 if (lp->scaling_used)
	  a/=(lp->scale[varnr]*lp->scale[0]);
	 if (prow[l]*(lp->lower[l] ? 1 : -1)<0.0) {
	  if (a<min1) min1=a;
	 }
	 else {
	  if (a<min2) min2=a;
	 }
	}
       if (!lp->lower[varnr]) {
	a=min1;
	min1=min2;
	min2=a;
       }
       if (min1<infinite) from=OrigObj[i]-min1;
       if (min2<infinite) till=OrigObj[i]+min2;
       if (lp->solution[varnr]==0.0) till=0.0; /* if value is 0 then there can't be an upper range */
      }
     }
     lp->objfrom[i]=from;
     lp->objtill[i]=till;
    }

    free(prow);
    free(OrigObj);
    free(drow);
  }
  return(ok);
 } /* calculate_sensitivity_obj */

static MYBOOL check_if_less(lprec *lp,
			  REAL x,
			  REAL y,
			  REAL value)
{
  if(x >= y) {
    report(lp, DETAILED, "Error: new upper or lower bound is not more restrictive");
    report(lp, DETAILED, "bound 1: %g, bound 2: %g, value: %g",
  	    (double)x, (double)y, (double)value);
    return(FALSE);
  }
  return(TRUE);
}

#if defined CHECK_SOLUTION
static short check_solution(lprec *lp,
                          int  lastcolumn,
		          REAL *solution,
			  REAL *upbo,
			  REAL *lowbo)
{
  REAL  test, value;
  int i, n;

  report(lp, NORMAL, "lp_solve successful at iteration %d with a best value of (%g)",
         lp->total_iter, solution[0]);
  if(lp->total_nodes > 1)
    report(lp, NORMAL, "lp_solve explored %d nodes to find optimum",
           lp->total_nodes);

/* Check if solution values are within the bounds; allowing a margin for numerical errors */
  n = 0;
  for(i = lp->rows + 1; i <= lp->rows+lastcolumn; i++) {
    test = lowbo[i];
    if(lp->columns_scaled)
      test *= lp->scale[i];

    if((solution[i] < test-SOLUTIONEPS*(1+fabs(test))) &&
	   !(lp->var_is_sc[i - lp->rows] > 0))  {
      report(lp, NORMAL,
              "Error: variable %s has a solution (%g) smaller than its lower bound (%g)",
              get_col_name(lp, i-lp->rows), (double)solution[i], (double)test);
      n++;
    }

    test = upbo[i];
    if(lp->columns_scaled)
      test *= lp->scale[i];

    if(solution[i] > test+SOLUTIONEPS*(1+fabs(test))) {
      report(lp, NORMAL,
              "Error: variable %s has a solution (%g) larger than its upper bound (%g)",
              get_col_name(lp, i-lp->rows), (double)solution[i], (double)test);
      n++;
    }
	if(n >= 20)
	  break;
  }

/* Check if constraint values are within the bounds; allowing a margin for numerical errors */
  for(i = 1; i <= lp->rows; i++) {

    value = solution[i];

    test = lp->orig_rh[i];
    if(lp->ch_sign[i]) {
      test = -test;
      test += fabs(upbo[i]);
    }
    if(lp->scaling_used)
      test /= lp->scale[i];

    if(value > test+SOLUTIONEPS*(1+fabs(test))) {
      report(lp, NORMAL,
              "Error: constraint %s has a value (%g) larger than its upper bound (%g)",
              get_row_name(lp, i), (double)value, (double)test);
      n++;
    }

    if(lp->ch_sign[i]) {
      test = lp->orig_rh[i];
	  test = -test;
    }
    else
      test = lp->orig_rh[i]-fabs(upbo[i]);
    if(lp->scaling_used)
      test /= lp->scale[i];

    if(value < test-SOLUTIONEPS*(1+fabs(test))) {
      report(lp, NORMAL,
              "Error: constraint %s has a value (%g) smaller than its lower bound (%g)",
              get_row_name(lp, i), (double)value, (double)test);
      n++;
    }

    if(n >= 20)
      break;
  }

  if(n > 0)
    return(FAILURE);
  else
    return(OPTIMAL);

} /* check_solution */
#endif /* CHECK_SOLUTION */

/* set working lower bounds to zero and transform rh correspondingly */
static void basetozero(lprec *lp)
{
  int i, j;
  LREAL theta;

  for(i = 1; i <= lp->columns; i++) {
    j = lp->rows + i;
    theta = lp->lowbo[j];
    if(theta != 0) {
      if(lp->upbo[j] < lp->infinite)
        lp->upbo[j] = (REAL) (lp->upbo[j] - theta);
      for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++)
        lp->rh[lp->mat[j].row_nr] = (REAL) (lp->rh[lp->mat[j].row_nr] - theta * lp->mat[j].value);
    }
  }
  /* Added by KE */
/*  for(i = 1; i <= lp->rows; i++)
    my_round(lp->rh[i], lp->epsel); */
}

static REAL int_floor(lprec *lp, int column, REAL value)
{
  value = floor(value);
  if(lp->columns_scaled && (lp->scalemode & INTEGERSCALE)) {
    value /= lp->scale[column];
    value += lp->epsel;
  }
  return(value);
}

static REAL int_ceil(lprec *lp, int column, REAL value)
{
  value = ceil(value);
  if(lp->columns_scaled && (lp->scalemode & INTEGERSCALE)) {
    value /= lp->scale[column];
    value -= lp->epsel;
  }
  return(value);
}

static short branch_and_bound(lprec *lp, REAL *upbo, REAL *lowbo, int notint, MYBOOL prunemode)
{
/* set up two new problems for the normal case.  If prune is non-negative, however,
   the brancing gets truncated.  The floor is skipped if prune == 0,
   and the ceiling is skipped if prune > 0 */

  REAL   *new_upbo, *new_lowbo;
  REAL    new_bound;
  REAL    tmpreal, sc_bound;
  MYBOOL   *new_lower, *new_basis, ActiveSOS, IsSOS, IntegerSOS;
  int    *new_bas;
  int     i, ii, k, MILPCount;
  short   spx_saved, failure = RUNNING, resone = RUNNING, restwo = RUNNING;

  spx_saved = lp->spx_status;
  lp->spx_status = RUNNING;
  if(yieldformessages(lp)!=0)
    lp->spx_status = USERABORT;

  if((lp->usermessage != NULL) && (lp->msgmask & MSG_MILPSTRATEGY))
    lp->usermessage(lp, lp->msghandle, MSG_MILPSTRATEGY);

  if(lp->spx_status != RUNNING)
    return(lp->spx_status);
  lp->spx_status = spx_saved;

 /* allocate room for them */
  new_upbo = new_lowbo = NULL;
  new_lower = new_basis = NULL;
  new_bas = NULL;
  if ((MALLOCCPY(new_upbo, upbo, lp->sum + 1) == NULL) ||
      (MALLOCCPY(new_lowbo, lowbo, lp->sum + 1) == NULL) ||
      (MALLOCCPY(new_lower, lp->lower, lp->sum + 1) == NULL) ||
      (MALLOCCPY(new_basis, lp->basis, lp->sum + 1) == NULL) ||
      (MALLOCCPY(new_bas, lp->bas, lp->rows + 1) == NULL)) {
    FREE(new_bas);
    FREE(new_basis);
    FREE(new_lower);
    FREE(new_lowbo);
    FREE(new_upbo);
    return(OUT_OF_MEMORY);
  }

 /* set local pruning info, automatic, or user-defined strategy */
  if(prunemode != AUTOMATIC) {
    i = prunemode;
  }
  else if(lp->floor_first == AUTOMATIC) {
    tmpreal = lp->solution[notint];
    tmpreal = tmpreal - (REAL)floor((double)tmpreal);
    if(tmpreal >= 0.5)
      i = FALSE;
    else
      i = TRUE;
    if((lp->bb_rule != BEST_SELECT) && (lp->bb_rule != GREEDY_SELECT))
      i = 1 - i;
  }
  else
    i = lp->floor_first;

 /* Force two runs through the MILP tree */
  resone = INFEASIBLE;
  restwo = INFEASIBLE;

  if(prunemode == TRUE)
    MILPCount = 1;
  else
    MILPCount = 2;

  tmpreal = lp->solution[notint];
  k = notint - lp->rows;
  IsSOS = (MYBOOL) (lp->sos_count > 0 && SOS_is_member(lp, 0, k));
  ActiveSOS = (MYBOOL) (IsSOS && prunemode == FALSE);
  IntegerSOS = (MYBOOL) (IsSOS && prunemode == AUTOMATIC);

  sc_bound = lp->var_is_sc[k];
  if(lp->scaling_used)
    sc_bound *= lp->scale[notint];

 /* Must make sure that we handle fractional lower bounds properly;
    also to ensure that we do a full binary tree search */
  if(is_int(lp,k) && ((sc_bound > 0) &&
                                      (tmpreal > floor(sc_bound)))) {
    tmpreal += 1;
    lowbo[notint] = int_floor(lp, notint, tmpreal);
  }

 /* SC logic: If the current SC variable value is in the [0..NZLOBOUND> range, then

	      UP: Set lower bound to NZLOBOUND, upper bound is the original
		  LO: Fix the variable by setting upper/lower bound to zero

    ... indicate that the variable is B&B-active by reversing sign of var_is_sc[]. */

  while (MILPCount) {
    if(i) {
      if((sc_bound > 0) && (tmpreal < sc_bound))
        new_bound = 0;
      else if(ActiveSOS) {
        new_bound = lp->orig_lowbo[notint];
      }
      else if(is_int(lp,k) && prunemode == AUTOMATIC)
        new_bound = int_floor(lp, notint, tmpreal);
      else if(lp->scaling_used)
  	new_bound = sc_bound/lp->scale[notint];
      else
  	new_bound = sc_bound;

    /* this bound might conflict */
      if(new_bound < lowbo[notint]) {
        debug_print(lp,
            "New upper bound value %g conflicts with old lower bound %g",
            (double)new_bound, (double)lowbo[notint]);
        resone = MILP_FAIL;
      }
      else { /* bound feasible */
        if(lp->debug)  /* Added by KE */
          check_if_less(lp, new_bound, upbo[notint], lp->solution[notint]);
        if(prunemode == TRUE) {
	  ii = 0+1;
	  if((SOS_fix_unmarked(lp, k, 0, new_upbo, new_bound, TRUE, &ii) < 0) ||
	     (ii == 0)) {
            MILPCount--;
            i = FALSE;
          }
        }
	else if(IsSOS && new_bound != 0 && !IntegerSOS) {
          MILPCount--;
          i = FALSE;
        }
        else
          new_upbo[notint] = new_bound;
        if (i) {
          debug_print(lp, "starting floor subproblem with bounds:");
          debug_print_bounds(lp, new_upbo, lowbo);

          if(IsSOS) SOS_set_marked(lp, 0, k, FALSE);
	  lp->var_is_sc[k] *= -1;
          lp->eta_valid = FALSE;

          resone = milpsolve(lp, new_upbo, lowbo, new_basis, new_lower, new_bas, TRUE);

          lp->eta_valid = FALSE;
          lp->var_is_sc[k] *= -1;
	  if(IsSOS) SOS_unmark(lp, 0, k, FALSE);
        }
      }
      if (i) {
        if((resone == USERABORT) || (resone == TIMEOUT) || (resone == OUT_OF_MEMORY)) {
          failure = resone;
          break;
        }
        MILPCount--;
        i = FALSE;
      }
    }
    else {
      if((sc_bound > 0) && (tmpreal < sc_bound)) {
        if(lp->scaling_used)
          new_bound = sc_bound / lp->scale[notint];
        else
	  new_bound = sc_bound;
	if(is_int(lp, k))
          new_bound = int_ceil(lp, notint, new_bound);
      }
      else if(ActiveSOS) {
        if(SOS_is_member_of_type(lp, k, SOS3))
          new_bound = 1;
	else
	  new_bound = lp->orig_lowbo[notint];
      }
      else if(is_int(lp,k) && prunemode == AUTOMATIC)
        new_bound = int_ceil(lp, notint, tmpreal);
      else {
	if(lp->scaling_used)
	  new_bound = sc_bound/lp->scale[notint];
        else
          new_bound = sc_bound;
      }

      if(new_bound > upbo[notint]) {
        debug_print(lp,
            "New lower bound value %g conflicts with old upper bound %g",
            (double)new_bound, (double)upbo[notint]);
        restwo = MILP_FAIL;
      }
      else { /* bound feasible */
        if(lp->debug)  /* Added by KE */
          check_if_less(lp, lowbo[notint], new_bound, lp->solution[notint]);

        new_lowbo[notint] = new_bound;
        debug_print(lp, "starting ceiling subproblem with bounds:");
        debug_print_bounds(lp, upbo, new_lowbo);

        if(IsSOS) SOS_set_marked(lp, 0, k, TRUE);
	lp->var_is_sc[k] *= -1;
        lp->eta_valid = FALSE;


/*	if(ActiveSOS) {
	  ii = SOS_fix_left(lp, k, 0, new_upbo, 0, TRUE);
	  if(ii < 0) {
            MILPCount--;
            i = FALSE;
            goto MILPRedo;
          }
          restwo = milpsolve(lp, new_upbo, new_lowbo, new_basis, new_lower, new_bas, TRUE);
        }
	else */
        restwo = milpsolve(lp, upbo, new_lowbo, new_basis, new_lower, new_bas, TRUE);

        lp->eta_valid = FALSE;
        lp->var_is_sc[k] *= -1;
        if(IsSOS) SOS_unmark(lp, 0, k, TRUE);
      }
      if((restwo == USERABORT) || (restwo == TIMEOUT) || (restwo == OUT_OF_MEMORY)) {
        failure = restwo;
        break;
      }
      MILPCount--;
      i = TRUE;
    }
  }

  if (MILPCount == 0) {
    /* Check for user abort or timout (resone was tested earlier) */
    if((resone != OPTIMAL) && (restwo != OPTIMAL)) /* both failed and must have been infeasible */
      failure = INFEASIBLE;
    else
      failure = OPTIMAL;
  }

  free(new_upbo);
  free(new_lowbo);
  free(new_basis);
  free(new_lower);
  free(new_bas);
  return(failure);
}

static short milpsolve(lprec *lp,
             REAL   *upbo,
             REAL   *lowbo,
             MYBOOL   *sbasis,
             MYBOOL   *slower,
             int    *sbas,
             int     recursive)
{
  int i, j, k, tilted, is_better, is_equal;
  int notint, nr_not_int;
  short failure;
  MYBOOL is_sos_int;
  REAL tmpreal;
  MYBOOL redo;

  if(lp->Break_bb)
    return(BREAK_BB);

  lp->Level++;
  lp->total_nodes++;

  if(lp->Level > lp->max_level)
    lp->max_level = lp->Level;

  debug_print(lp, "starting milpsolve");

  /* make fresh copies of upbo, lowbo, rh as solving changes them */
  /* Converted to macro by KE */
  MEMCPY(lp->upbo,  upbo,        lp->sum + 1);
  MEMCPY(lp->lowbo, lowbo,       lp->sum + 1);
  MEMCPY(lp->rh,    lp->orig_rh, lp->rows + 1);

  /* make sure we do not do memcpy(lp->basis, lp->basis ...) ! */
  /* Converted to macro by KE */
  if(recursive) {
    MEMCPY(lp->basis, sbasis, lp->sum + 1);
    MEMCPY(lp->lower, slower, lp->sum