Files
ortools-clone/ortools/glop/lp_solver.h
2024-01-04 13:43:15 +01:00

312 lines
14 KiB
C++

// Copyright 2010-2024 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef OR_TOOLS_GLOP_LP_SOLVER_H_
#define OR_TOOLS_GLOP_LP_SOLVER_H_
#include <memory>
#include <string>
#include "ortools/glop/parameters.pb.h"
#include "ortools/glop/revised_simplex.h"
#include "ortools/lp_data/lp_data.h"
#include "ortools/lp_data/lp_types.h"
#include "ortools/util/logging.h"
#include "ortools/util/time_limit.h"
namespace operations_research {
namespace glop {
// A full-fledged linear programming solver.
class LPSolver {
public:
LPSolver();
// This type is neither copyable nor movable.
LPSolver(const LPSolver&) = delete;
LPSolver& operator=(const LPSolver&) = delete;
// Sets and gets the solver parameters.
// See the proto for an extensive documentation.
void SetParameters(const GlopParameters& parameters);
const GlopParameters& GetParameters() const;
GlopParameters* GetMutableParameters();
// Returns a string that describes the version of the solver.
static std::string GlopVersion();
// Solves the given linear program and returns the solve status. See the
// ProblemStatus documentation for a description of the different values.
//
// The solution can be retrieved afterwards using the getter functions below.
// Note that depending on the returned ProblemStatus the solution values may
// not mean much, so it is important to check the returned status.
//
// Incrementality: From one Solve() call to the next, the internal state is
// not cleared and the solver may take advantage of its current state if the
// given lp is only slightly modified. If the modification is too important,
// or if the solver does not see how to reuse the previous state efficiently,
// it will just solve the problem from scratch. On the other hand, if the lp
// is the same, calling Solve() again should basically resume the solve from
// the last position. To disable this behavior, simply call Clear() before.
ABSL_MUST_USE_RESULT ProblemStatus Solve(const LinearProgram& lp);
// Same as Solve() but use the given time limit rather than constructing a new
// one from the current GlopParameters.
ABSL_MUST_USE_RESULT ProblemStatus SolveWithTimeLimit(const LinearProgram& lp,
TimeLimit* time_limit);
// Puts the solver in a clean state.
//
// Calling Solve() for the first time, or calling Clear() then Solve() on the
// same problem is guaranteed to be deterministic and to always give the same
// result, assuming that no time limit was specified.
void Clear();
// Advanced usage. This should be called before calling Solve(). It will
// configure the solver to try to start from the given point for the next
// Solve() only. Note that calling Clear() will invalidate this information.
//
// If the set of variables/constraints with a BASIC status does not form a
// basis a warning will be logged and the code will ignore it. Otherwise, the
// non-basic variables will be initialized to their given status and solving
// will start from there (even if the solution is not primal/dual feasible).
//
// Important: There is no facility to transform this information in sync with
// presolve. So you should probably disable presolve when using this since
// otherwise there is a good chance that the matrix will change and that the
// given basis will make no sense. Even worse if it happens to be factorizable
// but doesn't correspond to what was intended.
void SetInitialBasis(const VariableStatusRow& variable_statuses,
const ConstraintStatusColumn& constraint_statuses);
// This loads a given solution and computes related quantities so that the
// getters below will refer to it.
//
// Depending on the given solution status, this also checks the solution
// feasibility or optimality. The exact behavior and tolerances are controlled
// by the solver parameters. Because of this, the returned ProblemStatus may
// be changed from the one passed in the ProblemSolution to ABNORMAL or
// IMPRECISE. Note that this is the same logic as the one used by Solve() to
// verify the solver solution.
ProblemStatus LoadAndVerifySolution(const LinearProgram& lp,
const ProblemSolution& solution);
// Returns the objective value of the solution with its offset and scaling.
Fractional GetObjectiveValue() const;
// Accessors to information related to variables.
const DenseRow& variable_values() const { return primal_values_; }
const DenseRow& reduced_costs() const { return reduced_costs_; }
const VariableStatusRow& variable_statuses() const {
return variable_statuses_;
}
// Accessors to information related to constraints. The activity of a
// constraint is the sum of its linear terms evaluated with variables taking
// their values at the current solution.
//
// Note that the dual_values() do not take into account an eventual objective
// scaling of the solved LinearProgram.
const DenseColumn& dual_values() const { return dual_values_; }
const DenseColumn& constraint_activities() const {
return constraint_activities_;
}
const ConstraintStatusColumn& constraint_statuses() const {
return constraint_statuses_;
}
// Accessors to information related to unboundedness. A primal ray is returned
// for primal unbounded problems and a dual ray is returned for dual unbounded
// problems. constraints_dual_ray corresponds to dual multiplier for
// constraints and variable_bounds_dual_ray corresponds to dual multipliers
// for variable bounds (cf. reduced_costs).
const DenseRow& primal_ray() const { return primal_ray_; }
const DenseColumn& constraints_dual_ray() const {
return constraints_dual_ray_;
}
const DenseRow& variable_bounds_dual_ray() const {
return variable_bounds_dual_ray_;
}
// Returns the primal maximum infeasibility of the solution.
// This indicates by how much the variable and constraint bounds are violated.
Fractional GetMaximumPrimalInfeasibility() const;
// Returns the dual maximum infeasibility of the solution.
// This indicates by how much the variable costs (i.e. objective) should be
// modified for the solution to be an exact optimal solution.
Fractional GetMaximumDualInfeasibility() const;
// Returns true if the solution status was OPTIMAL and it seems that there is
// more than one basic optimal solution. Note that this solver always returns
// an optimal BASIC solution and that there is only a finite number of them.
// Moreover, given one basic solution, since the basis is always refactorized
// at optimality before reporting the numerical result, then all the
// quantities (even the floating point ones) should be always the same.
//
// TODO(user): Test this behavior extensively if a client relies on it.
bool MayHaveMultipleOptimalSolutions() const;
// Returns the number of simplex iterations used by the last Solve().
int GetNumberOfSimplexIterations() const;
// Returns the "deterministic time" since the creation of the solver. Note
// That this time is only increased when some operations take place in this
// class.
//
// TODO(user): Currently, this is only modified when the simplex code is
// executed.
//
// TODO(user): Improve the correlation with the running time.
double DeterministicTime() const;
// Returns the SolverLogger used during solves.
//
// Please note that EnableLogging() and SetLogToStdOut() are reset at the
// beginning of each solve based on parameters so setting them will have no
// effect.
SolverLogger& GetSolverLogger();
private:
// Resizes all the solution vectors to the given sizes.
// This is used in case of error to make sure all the getter functions will
// not crash when given row/col inside the initial linear program dimension.
void ResizeSolution(RowIndex num_rows, ColIndex num_cols);
// Make sure the primal and dual values are within their bounds in order to
// have a strong guarantee on the optimal solution. See
// provide_strong_optimal_guarantee in the GlopParameters proto.
void MovePrimalValuesWithinBounds(const LinearProgram& lp);
void MoveDualValuesWithinBounds(const LinearProgram& lp);
// Runs the revised simplex algorithm if needed (i.e. if the program was not
// already solved by the preprocessors).
void RunRevisedSimplexIfNeeded(ProblemSolution* solution,
TimeLimit* time_limit);
// Checks that the returned solution values and statuses are consistent.
// Returns true if this is the case. See the code for the exact check
// performed.
bool IsProblemSolutionConsistent(const LinearProgram& lp,
const ProblemSolution& solution) const;
// Returns true if there may be multiple optimal solutions.
// The return value is true if:
// - a non-fixed variable, at one of its boumds, has its reduced
// cost close to zero.
// or if:
// - a non-equality constraint (i.e. l <= a.x <= r, with l != r), is at one of
// its bounds (a.x = r or a.x = l) and has its dual value close to zero.
bool IsOptimalSolutionOnFacet(const LinearProgram& lp);
// Computes derived quantities from the solution.
void ComputeReducedCosts(const LinearProgram& lp);
void ComputeConstraintActivities(const LinearProgram& lp);
// Computes the primal/dual objectives (without the offset). Note that the
// dual objective needs the reduced costs in addition to the dual values.
double ComputeObjective(const LinearProgram& lp);
double ComputeDualObjective(const LinearProgram& lp);
// Given a relative precision on the primal values of up to
// solution_feasibility_tolerance(), this returns an upper bound on the
// expected precision of the objective.
double ComputeMaxExpectedObjectiveError(const LinearProgram& lp);
// Returns the max absolute cost pertubation (resp. rhs perturbation) so that
// the pair (primal values, dual values) is an EXACT optimal solution to the
// perturbed problem. Note that this assumes that
// MovePrimalValuesWithinBounds() and MoveDualValuesWithinBounds() have
// already been called. The Boolean is_too_large is set to true if any of the
// perturbation exceed the tolerance (which depends of the coordinate).
//
// These bounds are computed using the variable and constraint statuses by
// enforcing the complementary slackness optimal conditions. Note that they
// are almost the same as ComputeActivityInfeasibility() and
// ComputeReducedCostInfeasibility() but looks for optimality rather than just
// feasibility.
//
// Note(user): We could get EXACT bounds on these perturbations by changing
// the rounding mode appropriately during these computations. But this is
// probably not needed.
Fractional ComputeMaxCostPerturbationToEnforceOptimality(
const LinearProgram& lp, bool* is_too_large);
Fractional ComputeMaxRhsPerturbationToEnforceOptimality(
const LinearProgram& lp, bool* is_too_large);
// Computes the maximum of the infeasibilities associated with each values.
// The returned infeasibilities are the maximum of the "absolute" errors of
// each vector coefficients.
//
// These function also set is_too_large to true if any infeasibility is
// greater than the tolerance (which depends of the coordinate).
double ComputePrimalValueInfeasibility(const LinearProgram& lp,
bool* is_too_large);
double ComputeActivityInfeasibility(const LinearProgram& lp,
bool* is_too_large);
double ComputeDualValueInfeasibility(const LinearProgram& lp,
bool* is_too_large);
double ComputeReducedCostInfeasibility(const LinearProgram& lp,
bool* is_too_large);
// On a call to Solve(), this is initialized to an exact copy of the given
// linear program. It is later modified by the preprocessors and then solved
// by the revised simplex.
//
// This is not efficient memory-wise but allows to check optimality with
// respect to the given LinearProgram that is guaranteed to not have been
// modified. It also allows for a nicer Solve() API with a const
// LinearProgram& input.
LinearProgram current_linear_program_;
SolverLogger logger_;
// The revised simplex solver.
std::unique_ptr<RevisedSimplex> revised_simplex_;
// The number of revised simplex iterations used by the last Solve().
int num_revised_simplex_iterations_;
// The current ProblemSolution.
// TODO(user): use a ProblemSolution directly? Note, that primal_ray_,
// constraints_dual_ray_ and variable_bounds_dual_ray_ are not currently in
// ProblemSolution and are filled directly by RunRevisedSimplexIfNeeded().
DenseRow primal_values_;
DenseColumn dual_values_;
VariableStatusRow variable_statuses_;
ConstraintStatusColumn constraint_statuses_;
DenseRow primal_ray_;
DenseColumn constraints_dual_ray_;
DenseRow variable_bounds_dual_ray_;
// Quantities computed from the solution and the linear program.
DenseRow reduced_costs_;
DenseColumn constraint_activities_;
Fractional problem_objective_value_ = 0.0;
bool may_have_multiple_solutions_;
Fractional max_absolute_primal_infeasibility_;
Fractional max_absolute_dual_infeasibility_;
// Proto holding all the parameters of the algorithm.
GlopParameters parameters_;
// The number of times Solve() was called. Used to number dump files.
int num_solves_;
};
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_LP_SOLVER_H_