forked from google/or-tools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
lp_solver.h
309 lines (266 loc) · 13.8 KB
/
lp_solver.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
// Copyright 2010-2022 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();
// 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_;
DISALLOW_COPY_AND_ASSIGN(LPSolver);
};
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_LP_SOLVER_H_