Google OR-Tools v9.12
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
revised_simplex.h
Go to the documentation of this file.
1// Copyright 2010-2025 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
14// Solves a Linear Programming problem using the Revised Simplex algorithm
15// as described by G.B. Dantzig.
16// The general form is:
17// min c.x where c and x are n-vectors,
18// subject to Ax = b where A is an mxn-matrix, b an m-vector,
19// with l <= x <= u, i.e.
20// l_i <= x_i <= u_i for all i in {1 .. m}.
21//
22// c.x is called the objective function.
23// Each row a_i of A is an n-vector, and a_i.x = b_i is a linear constraint.
24// A is called the constraint matrix.
25// b is called the right hand side (rhs) of the problem.
26// The constraints l_i <= x_i <= u_i are called the generalized bounds
27// of the problem (most introductory textbooks only deal with x_i >= 0, as
28// did the first version of the Simplex algorithm). Note that l_i and u_i
29// can be -infinity and +infinity, respectively.
30//
31// To simplify the entry of data, this code actually handles problems in the
32// form:
33// min c.x where c and x are n-vectors,
34// subject to:
35// A1 x <= b1
36// A2 x >= b2
37// A3 x = b3
38// l <= x <= u
39//
40// It transforms the above problem into
41// min c.x where c and x are n-vectors,
42// subject to:
43// A1 x + s1 = b1
44// A2 x - s2 = b2
45// A3 x = b3
46// l <= x <= u
47// s1 >= 0, s2 >= 0
48// where xT = (x1, x2, x3),
49// s1 is an m1-vector (m1 being the height of A1),
50// s2 is an m2-vector (m2 being the height of A2).
51//
52// The following are very good references for terminology, data structures,
53// and algorithms. They all contain a wealth of references.
54//
55// Vasek Chvátal, "Linear Programming," W.H. Freeman, 1983. ISBN 978-0716715870.
56// http://www.amazon.com/dp/0716715872
57//
58// Robert J. Vanderbei, "Linear Programming: Foundations and Extensions,"
59// Springer, 2010, ISBN-13: 978-1441944979
60// http://www.amazon.com/dp/1441944974
61//
62// Istvan Maros, "Computational Techniques of the Simplex Method.", Springer,
63// 2002, ISBN 978-1402073328
64// http://www.amazon.com/dp/1402073321
65//
66// ===============================================
67// Short description of the dual simplex algorithm.
68//
69// The dual simplex algorithm uses the same data structure as the primal, but
70// progresses towards the optimal solution in a different way:
71// * It tries to keep the dual values dual-feasible at all time which means that
72// the reduced costs are of the correct sign depending on the bounds of the
73// non-basic variables. As a consequence the values of the basic variable are
74// out of bound until the optimal is reached.
75// * A basic leaving variable is selected first (dual pricing) and then a
76// corresponding entering variable is selected. This is done in such a way
77// that the dual objective value increases (lower bound on the optimal
78// solution).
79// * Once the basis pivot is chosen, the variable values and the reduced costs
80// are updated the same way as in the primal algorithm.
81//
82// Good references on the Dual simplex algorithm are:
83//
84// Robert Fourer, "Notes on the Dual simplex Method", March 14, 1994.
85// http://users.iems.northwestern.edu/~4er/WRITINGS/dual.pdf
86//
87// Achim Koberstein, "The dual simplex method, techniques for a fast and stable
88// implementation", PhD, Paderborn, Univ., 2005.
89// http://digital.ub.uni-paderborn.de/hs/download/pdf/3885?originalFilename=true
90
91#ifndef OR_TOOLS_GLOP_REVISED_SIMPLEX_H_
92#define OR_TOOLS_GLOP_REVISED_SIMPLEX_H_
93
94#include <cstdint>
95#include <string>
96#include <vector>
97
98#include "absl/base/attributes.h"
99#include "absl/log/die_if_null.h"
100#include "absl/random/bit_gen_ref.h"
101#include "absl/random/random.h"
102#include "ortools/base/types.h"
107#include "ortools/glop/parameters.pb.h"
108#include "ortools/glop/pricing.h"
111#include "ortools/glop/status.h"
122#include "ortools/util/logging.h"
124#include "ortools/util/stats.h"
126
127namespace operations_research {
128namespace glop {
129
130// Entry point of the revised simplex algorithm implementation.
131class RevisedSimplex {
132 public:
134
135 // This type is neither copyable nor movable.
136 RevisedSimplex(const RevisedSimplex&) = delete;
139 // Sets or gets the algorithm parameters to be used on the next Solve().
140 void SetParameters(const GlopParameters& parameters);
141 const GlopParameters& GetParameters() const { return parameters_; }
143 // Solves the given linear program.
144 //
145 // We accept two forms of LinearProgram:
146 // - The lp can be in the equations form Ax = 0 created by
147 // LinearProgram::AddSlackVariablesForAllRows(), i.e. the rightmost square
148 // submatrix of A is an identity matrix, all its columns have been marked as
149 // slack variables, and the bounds of all constraints have been set to 0.
150 // - If not, we will convert it internally while copying it to the internal
151 // structure used.
152 //
153 // By default, the algorithm tries to exploit the computation done during the
154 // last Solve() call. It will analyze the difference of the new linear program
155 // and try to use the previously computed solution as a warm-start. To disable
156 // this behavior or give explicit warm-start data, use one of the State*()
157 // functions below.
158 ABSL_MUST_USE_RESULT Status Solve(const LinearProgram& lp,
160
161 // Do not use the current solution as a warm-start for the next Solve(). The
162 // next Solve() will behave as if the class just got created.
164
165 // Uses the given state as a warm-start for the next Solve() call.
166 void LoadStateForNextSolve(const BasisState& state);
167
168 // Advanced usage. While constructing the initial basis, if this is called
169 // then we will use these values as the initial starting value for the FREE
170 // variables.
172
173 // Getters to retrieve all the information computed by the last Solve().
174 RowIndex GetProblemNumRows() const;
175 ColIndex GetProblemNumCols() const;
178 int64_t GetNumberOfIterations() const;
179 Fractional GetVariableValue(ColIndex col) const;
180 Fractional GetReducedCost(ColIndex col) const;
181 const DenseRow& GetReducedCosts() const;
182 Fractional GetDualValue(RowIndex row) const;
183 Fractional GetConstraintActivity(RowIndex row) const;
184 VariableStatus GetVariableStatus(ColIndex col) const;
185 ConstraintStatus GetConstraintStatus(RowIndex row) const;
186 const BasisState& GetState() const;
187 double DeterministicTime() const;
188 bool objective_limit_reached() const { return objective_limit_reached_; }
191 return dual_edge_norms_.GetEdgeSquaredNorms();
192 }
193
194 const DenseBitRow& GetNotBasicBitRow() const {
195 return variables_info_.GetNotBasicBitRow();
196 }
197
198 // If the problem status is PRIMAL_UNBOUNDED (respectively DUAL_UNBOUNDED),
199 // then the solver has a corresponding primal (respectively dual) ray to show
200 // the unboundness. From a primal (respectively dual) feasible solution any
201 // positive multiple of this ray can be added to the solution and keep it
202 // feasible. Moreover, by doing so, the objective of the problem will improve
203 // and its magnitude will go to infinity.
204 //
205 // Note that when the problem is DUAL_UNBOUNDED, the dual ray is also known as
206 // the Farkas proof of infeasibility of the problem.
207 const DenseRow& GetPrimalRay() const;
208 const DenseColumn& GetDualRay() const;
209
210 // This is the "dual ray" linear combination of the matrix rows.
211 const DenseRow& GetDualRayRowCombination() const;
212
213 // Returns the index of the column in the basis and the basis factorization.
214 // Note that the order of the column in the basis is important since it is the
215 // one used by the various solve functions provided by the BasisFactorization
216 // class.
217 ColIndex GetBasis(RowIndex row) const;
218
219 const ScatteredRow& GetUnitRowLeftInverse(RowIndex row) {
220 return update_row_.ComputeAndGetUnitRowLeftInverse(row);
221 }
222
223 // Returns a copy of basis_ vector for outside applications (like cuts) to
224 // have the correspondence between rows and columns of the dictionary.
225 RowToColMapping GetBasisVector() const { return basis_; }
228
229 // Returns statistics about this class as a string.
230 std::string StatString();
231
232 // Computes the dictionary B^-1*N on-the-fly row by row. Returns the resulting
233 // matrix as a vector of sparse rows so that it is easy to use it on the left
234 // side in the matrix multiplication. Runs in O(num_non_zeros_in_matrix).
235 // TODO(user): Use row scales as well.
236 RowMajorSparseMatrix ComputeDictionary(const DenseRow* column_scales);
237
238 // Initializes the matrix for the given 'linear_program' and 'state' and
239 // computes the variable values for basic variables using non-basic variables.
240 void ComputeBasicVariablesForState(const LinearProgram& linear_program,
241 const BasisState& state);
242
243 // This is used in a MIP context to polish the final basis. We assume that the
244 // columns for which SetIntegralityScale() has been called correspond to
245 // integral variable once multiplied by the given factor.
246 void ClearIntegralityScales() { integrality_scale_.clear(); }
247 void SetIntegralityScale(ColIndex col, Fractional scale);
248
249 void SetLogger(SolverLogger* logger) { logger_ = logger; }
251 // Advanced usage. For fast incremental call to the solver, it is better not
252 // to use LinearProgram at all. This api allows to directly modify the
253 // internal data of glop and then call solve.
254 const CompactSparseMatrix& MatrixWithSlack() const { return compact_matrix_; }
256 transpose_was_changed_ = true;
257 return &transposed_matrix_;
258 }
260 return variables_info_.MutableLowerBounds();
261 }
263 return variables_info_.MutableUpperBounds();
264 }
266 const DenseRow& objective, Fractional objective_scaling_factor,
267 Fractional objective_offset, TimeLimit* time_limit);
268
269 private:
270 struct IterationStats : public StatsGroup {
271 IterationStats()
272 : StatsGroup("IterationStats"),
273 total("total", this),
274 normal("normal", this),
275 bound_flip("bound_flip", this),
276 refactorize("refactorize", this),
277 degenerate("degenerate", this),
278 num_dual_flips("num_dual_flips", this),
279 degenerate_run_size("degenerate_run_size", this) {}
280 TimeDistribution total;
281 TimeDistribution normal;
282 TimeDistribution bound_flip;
283 TimeDistribution refactorize;
284 TimeDistribution degenerate;
285 IntegerDistribution num_dual_flips;
286 IntegerDistribution degenerate_run_size;
287 };
288
289 struct RatioTestStats : public StatsGroup {
290 RatioTestStats()
291 : StatsGroup("RatioTestStats"),
292 bound_shift("bound_shift", this),
293 abs_used_pivot("abs_used_pivot", this),
294 abs_tested_pivot("abs_tested_pivot", this),
295 abs_skipped_pivot("abs_skipped_pivot", this),
296 direction_density("direction_density", this),
297 leaving_choices("leaving_choices", this),
298 num_perfect_ties("num_perfect_ties", this) {}
299 DoubleDistribution bound_shift;
300 DoubleDistribution abs_used_pivot;
301 DoubleDistribution abs_tested_pivot;
302 DoubleDistribution abs_skipped_pivot;
303 RatioDistribution direction_density;
304 IntegerDistribution leaving_choices;
305 IntegerDistribution num_perfect_ties;
306 };
307
308 enum class Phase { FEASIBILITY, OPTIMIZATION, PUSH };
309
310 enum class RefactorizationReason {
311 DEFAULT,
312 SMALL_PIVOT,
313 IMPRECISE_PIVOT,
314 NORM,
315 RC,
316 VAR_VALUES,
317 FINAL_CHECK
318 };
319
320 ABSL_MUST_USE_RESULT Status SolveInternal(double start_time, bool maximize,
321 const DenseRow& objective,
322 TimeLimit* time_limit);
323
324 // Propagates parameters_ to all the other classes that need it.
325 //
326 // TODO(user): Maybe a better design is for them to have a reference to a
327 // unique parameters object? It will clutter a bit more these classes'
328 // constructor though.
329 void PropagateParameters();
330
331 // Returns a string containing the same information as with GetSolverStats,
332 // but in a much more human-readable format. For example:
333 // Problem status : Optimal
334 // Solving time : 1.843
335 // Number of iterations : 12345
336 // Time for solvability (first phase) : 1.343
337 // Number of iterations for solvability : 10000
338 // Time for optimization : 0.5
339 // Number of iterations for optimization : 2345
340 // Maximum time allowed in seconds : 6000
341 // Maximum number of iterations : 1000000
342 // Stop after first basis : 0
343 std::string GetPrettySolverStats() const;
344
345 // Returns a string containing formatted information about the variable
346 // corresponding to column col.
347 std::string SimpleVariableInfo(ColIndex col) const;
348
349 // Displays a short string with the current iteration and objective value.
350 void DisplayIterationInfo(bool primal, RefactorizationReason reason =
351 RefactorizationReason::DEFAULT);
352
353 // Displays the error bounds of the current solution.
354 void DisplayErrors();
355
356 // Displays the status of the variables.
357 void DisplayInfoOnVariables() const;
358
359 // Displays the bounds of the variables.
360 void DisplayVariableBounds();
361
362 // Displays the following information:
363 // * Linear Programming problem as a dictionary, taking into
364 // account the iterations that have been made;
365 // * Variable info;
366 // * Reduced costs;
367 // * Variable bounds.
368 // A dictionary is in the form:
369 // xB = value + sum_{j in N} pa_ij x_j
370 // z = objective_value + sum_{i in N} rc_i x_i
371 // where the pa's are the coefficients of the matrix after the pivotings
372 // and the rc's are the reduced costs, i.e. the coefficients of the objective
373 // after the pivotings.
374 // Dictionaries are the modern way of presenting the result of an iteration
375 // of the Simplex algorithm in the literature.
376 void DisplayRevisedSimplexDebugInfo();
377
378 // Displays the Linear Programming problem as it was input.
379 void DisplayProblem() const;
380
381 // Returns the current objective value. This is just the sum of the current
382 // variable values times their current cost.
383 Fractional ComputeObjectiveValue() const;
384
385 // Returns the current objective of the linear program given to Solve() using
386 // the initial costs, maximization direction, objective offset and objective
387 // scaling factor.
388 Fractional ComputeInitialProblemObjectiveValue() const;
389
390 // Assigns names to variables. Variables in the input will be named
391 // x1..., slack variables will be s1... .
392 void SetVariableNames();
393
394 // Sets the variable status and derives the variable value according to the
395 // exact status definition. This can only be called for non-basic variables
396 // because the value of a basic variable is computed from the values of the
397 // non-basic variables.
398 void SetNonBasicVariableStatusAndDeriveValue(ColIndex col,
399 VariableStatus status);
400
401 // Checks if the basis_ and is_basic_ arrays are well formed. Also checks that
402 // the variable statuses are consistent with this basis. Returns true if this
403 // is the case. This is meant to be used in debug mode only.
404 bool BasisIsConsistent() const;
405
406 // Moves the column entering_col into the basis at position basis_row. Removes
407 // the current basis column at position basis_row from the basis and sets its
408 // status to leaving_variable_status.
409 void UpdateBasis(ColIndex entering_col, RowIndex basis_row,
410 VariableStatus leaving_variable_status);
411
412 // Initializes matrix-related internal data. Returns true if this data was
413 // unchanged. If not, also sets only_change_is_new_rows to true if compared
414 // to the current matrix, the only difference is that new rows have been
415 // added (with their corresponding extra slack variables). Similarly, sets
416 // only_change_is_new_cols to true if the only difference is that new columns
417 // have been added, in which case also sets num_new_cols to the number of
418 // new columns.
419 bool InitializeMatrixAndTestIfUnchanged(const LinearProgram& lp,
420 bool lp_is_in_equation_form,
421 bool* only_change_is_new_rows,
422 bool* only_change_is_new_cols,
423 ColIndex* num_new_cols);
424
425 // Checks if the only change to the bounds is the addition of new columns,
426 // and that the new columns have at least one bound equal to zero.
427 bool OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero(
428 const LinearProgram& lp, bool lp_is_in_equation_form,
429 ColIndex num_new_cols);
430
431 // Initializes objective-related internal data. Returns true if unchanged.
432 bool InitializeObjectiveAndTestIfUnchanged(const LinearProgram& lp);
433
434 // Computes the stopping criterion on the problem objective value.
435 void InitializeObjectiveLimit();
436
437 // Initializes the starting basis. In most cases it starts by the all slack
438 // basis and tries to apply some heuristics to replace fixed variables.
439 ABSL_MUST_USE_RESULT Status CreateInitialBasis();
440
441 // Sets the initial basis to the given columns, try to factorize it and
442 // recompute the basic variable values.
443 ABSL_MUST_USE_RESULT Status
444 InitializeFirstBasis(const RowToColMapping& initial_basis);
445
446 // Entry point for the solver initialization.
447 ABSL_MUST_USE_RESULT Status Initialize(const LinearProgram& lp);
448 ABSL_MUST_USE_RESULT Status FinishInitialization(bool solve_from_scratch);
449
450 // Saves the current variable statuses in solution_state_.
451 void SaveState();
452
453 // Displays statistics on what kinds of variables are in the current basis.
454 void DisplayBasicVariableStatistics();
455
456 // Tries to reduce the initial infeasibility (stored in error_) by using the
457 // singleton columns present in the problem. A singleton column is a column
458 // with only one non-zero. This is used by CreateInitialBasis().
459 void UseSingletonColumnInInitialBasis(RowToColMapping* basis);
460
461 // Returns the number of empty rows in the matrix, i.e. rows where all
462 // the coefficients are zero.
463 RowIndex ComputeNumberOfEmptyRows();
464
465 // Returns the number of empty columns in the matrix, i.e. columns where all
466 // the coefficients are zero.
467 ColIndex ComputeNumberOfEmptyColumns();
468
469 // Returns the number of super-basic variables. These are non-basic variables
470 // that are not at their bounds (if they have bounds), or non-basic free
471 // variables that are not at zero.
472 int ComputeNumberOfSuperBasicVariables() const;
473
474 // This method transforms a basis for the first phase, with the optimal
475 // value at zero, into a feasible basis for the initial problem, thus
476 // preparing the execution of phase-II of the algorithm.
477 void CleanUpBasis();
478
479 // If the primal maximum residual is too large, recomputes the basic variable
480 // value from the non-basic ones. This function also perturbs the bounds
481 // during the primal simplex if too many iterations are degenerate.
482 //
483 // Only call this on a refactorized basis to have the best precision.
484 void CorrectErrorsOnVariableValues();
485
486 // Computes b - A.x in error_
487 void ComputeVariableValuesError();
488
489 // Solves the system B.d = a where a is the entering column (given by col).
490 // Known as FTRAN (Forward transformation) in FORTRAN codes.
491 // See Chvatal's book for more detail (Chapter 7).
492 void ComputeDirection(ColIndex col);
493
494 // Computes a - B.d in error_ and return the maximum std::abs() of its coeffs.
495 Fractional ComputeDirectionError(ColIndex col);
496
497 // Computes the ratio of the basic variable corresponding to 'row'. A target
498 // bound (upper or lower) is chosen depending on the sign of the entering
499 // reduced cost and the sign of the direction 'd_[row]'. The ratio is such
500 // that adding 'ratio * d_[row]' to the variable value changes it to its
501 // target bound.
502 template <bool is_entering_reduced_cost_positive>
503 Fractional GetRatio(const DenseRow& lower_bounds,
504 const DenseRow& upper_bounds, RowIndex row) const;
505
506 // First pass of the Harris ratio test. Returns the harris ratio value which
507 // is an upper bound on the ratio value that the leaving variable can take.
508 // Fills leaving_candidates with the ratio and row index of a super-set of the
509 // columns with a ratio <= harris_ratio.
510 template <bool is_entering_reduced_cost_positive>
511 Fractional ComputeHarrisRatioAndLeavingCandidates(
512 Fractional bound_flip_ratio, SparseColumn* leaving_candidates) const;
513
514 // Chooses the leaving variable, considering the entering column and its
515 // associated reduced cost. If there was a precision issue and the basis is
516 // not refactorized, set refactorize to true. Otherwise, the row number of the
517 // leaving variable is written in *leaving_row, and the step length
518 // is written in *step_length.
519 Status ChooseLeavingVariableRow(ColIndex entering_col,
520 Fractional reduced_cost, bool* refactorize,
521 RowIndex* leaving_row,
522 Fractional* step_length,
523 Fractional* target_bound);
524
525 // Chooses the leaving variable for the primal phase-I algorithm. The
526 // algorithm follows more or less what is described in Istvan Maros's book in
527 // chapter 9.6 and what is done for the dual phase-I algorithm which was
528 // derived from Koberstein's PhD. Both references can be found at the top of
529 // this file.
530 void PrimalPhaseIChooseLeavingVariableRow(ColIndex entering_col,
531 Fractional reduced_cost,
532 bool* refactorize,
533 RowIndex* leaving_row,
534 Fractional* step_length,
535 Fractional* target_bound) const;
536
537 // Chooses an infeasible basic variable. The returned values are:
538 // - leaving_row: the basic index of the infeasible leaving variable
539 // or kNoLeavingVariable if no such row exists: the dual simplex algorithm
540 // has terminated and the optimal has been reached.
541 // - cost_variation: how much do we improve the objective by moving one unit
542 // along this dual edge.
543 // - target_bound: the bound at which the leaving variable should go when
544 // leaving the basis.
545 ABSL_MUST_USE_RESULT Status DualChooseLeavingVariableRow(
546 RowIndex* leaving_row, Fractional* cost_variation,
547 Fractional* target_bound);
548
549 // Updates the prices used by DualChooseLeavingVariableRow() after a simplex
550 // iteration by using direction_. The prices are stored in
551 // dual_pricing_vector_. Note that this function only takes care of the
552 // entering and leaving column dual feasibility status change and that other
553 // changes will be dealt with by DualPhaseIUpdatePriceOnReducedCostsChange().
554 void DualPhaseIUpdatePrice(RowIndex leaving_row, ColIndex entering_col);
555
556 // This must be called each time the dual_pricing_vector_ is changed at
557 // position row.
558 template <bool use_dense_update = false>
559 void OnDualPriceChange(DenseColumn::ConstView squared_norms, RowIndex row,
560 VariableType type, Fractional threshold);
561
562 // Updates the prices used by DualChooseLeavingVariableRow() when the reduced
563 // costs of the given columns have changed.
564 template <typename Cols>
565 void DualPhaseIUpdatePriceOnReducedCostChange(const Cols& cols);
566
567 // Same as DualChooseLeavingVariableRow() but for the phase I of the dual
568 // simplex. Here the objective is not to minimize the primal infeasibility,
569 // but the dual one, so the variable is not chosen in the same way. See
570 // "Notes on the Dual simplex Method" or Istvan Maros, "A Piecewise Linear
571 // Dual Phase-1 Algorithm for the Simplex Method", Computational Optimization
572 // and Applications, October 2003, Volume 26, Issue 1, pp 63-81.
573 // http://rd.springer.com/article/10.1023%2FA%3A1025102305440
574 ABSL_MUST_USE_RESULT Status DualPhaseIChooseLeavingVariableRow(
575 RowIndex* leaving_row, Fractional* cost_variation,
576 Fractional* target_bound);
577
578 // Makes sure the boxed variable are dual-feasible by setting them to the
579 // correct bound according to their reduced costs. This is called
580 // Dual feasibility correction in the literature.
581 //
582 // Note that this function is also used as a part of the bound flipping ratio
583 // test by flipping the boxed dual-infeasible variables at each iteration.
584 //
585 // If update_basic_values is true, the basic variable values are updated.
586 template <typename BoxedVariableCols>
587 void MakeBoxedVariableDualFeasible(const BoxedVariableCols& cols,
588 bool update_basic_values);
589
590 // Computes the step needed to move the leaving_row basic variable to the
591 // given target bound.
592 Fractional ComputeStepToMoveBasicVariableToBound(RowIndex leaving_row,
593 Fractional target_bound);
594
595 // Returns true if the basis obtained after the given pivot can be factorized.
596 bool TestPivot(ColIndex entering_col, RowIndex leaving_row);
597
598 // Gets the current LU column permutation from basis_representation,
599 // applies it to basis_ and then sets it to the identity permutation since
600 // it will no longer be needed during solves. This function also updates all
601 // the data that depends on the column order in basis_.
602 void PermuteBasis();
603
604 // Updates the system state according to the given basis pivot.
605 // Returns an error if the update could not be done because of some precision
606 // issue.
607 ABSL_MUST_USE_RESULT Status UpdateAndPivot(ColIndex entering_col,
608 RowIndex leaving_row,
609 Fractional target_bound);
610
611 // Displays all the timing stats related to the calling object.
612 void DisplayAllStats();
613
614 // Calls basis_factorization_.Refactorize() if refactorize is true, and
615 // returns its status. This also sets refactorize to false and invalidates any
616 // data structure that depends on the current factorization.
617 //
618 // The general idea is that if a refactorization is going to be needed during
619 // a simplex iteration, it is better to do it as soon as possible so that
620 // every component can take advantage of it.
621 Status RefactorizeBasisIfNeeded(bool* refactorize);
622
623 // Main iteration loop of the primal simplex.
624 ABSL_MUST_USE_RESULT Status PrimalMinimize(TimeLimit* time_limit);
625
626 // Main iteration loop of the dual simplex.
627 ABSL_MUST_USE_RESULT Status DualMinimize(bool feasibility_phase,
628 TimeLimit* time_limit);
629
630 // Pushes all super-basic variables to bounds (if applicable) or to zero (if
631 // unconstrained). This is part of a "crossover" procedure to find a vertex
632 // solution given a (near) optimal solution. Assumes that Minimize() or
633 // DualMinimize() has already run, i.e., that we are at an optimal solution
634 // within numerical tolerances.
635 ABSL_MUST_USE_RESULT Status PrimalPush(TimeLimit* time_limit);
636
637 // Experimental. This is useful in a MIP context. It performs a few degenerate
638 // pivot to try to mimize the fractionality of the optimal basis.
639 //
640 // We assume that the columns for which SetIntegralityScale() has been called
641 // correspond to integral variable once scaled by the given factor.
642 //
643 // I could only find slides for the reference of this "LP Solution Polishing
644 // to improve MIP Performance", Matthias Miltenberger, Zuse Institute Berlin.
645 ABSL_MUST_USE_RESULT Status Polish(TimeLimit* time_limit);
646
647 // Utility functions to return the current ColIndex of the slack column with
648 // given number. Note that currently, such columns are always present in the
649 // internal representation of a linear program.
650 ColIndex SlackColIndex(RowIndex row) const;
651
652 // Advances the deterministic time in time_limit with the difference between
653 // the current internal deterministic time and the internal deterministic time
654 // during the last call to this method.
655 // TODO(user): Update the internals of revised simplex so that the time
656 // limit is updated at the source and remove this method.
657 void AdvanceDeterministicTime(TimeLimit* time_limit);
658
659 // Problem status
660 ProblemStatus problem_status_;
661
662 // Current number of rows in the problem.
663 RowIndex num_rows_ = RowIndex(0);
664
665 // Current number of columns in the problem.
666 ColIndex num_cols_ = ColIndex(0);
667
668 // Index of the first slack variable in the input problem. We assume that all
669 // variables with index greater or equal to first_slack_col_ are slack
670 // variables.
671 ColIndex first_slack_col_ = ColIndex(0);
672
673 // We're using vectors after profiling and looking at the generated assembly
674 // it's as fast as std::unique_ptr as long as the size is properly reserved
675 // beforehand.
676
677 // Compact version of the matrix given to Solve().
678 CompactSparseMatrix compact_matrix_;
679
680 // The transpose of compact_matrix_, it may be empty if it is not needed.
681 CompactSparseMatrix transposed_matrix_;
682
683 // Stop the algorithm and report feasibility if:
684 // - The primal simplex is used, the problem is primal-feasible and the
685 // current objective value is strictly lower than primal_objective_limit_.
686 // - The dual simplex is used, the problem is dual-feasible and the current
687 // objective value is strictly greater than dual_objective_limit_.
688 Fractional primal_objective_limit_;
689 Fractional dual_objective_limit_;
690
691 // Current objective (feasibility for Phase-I, user-provided for Phase-II).
692 DenseRow current_objective_;
693
694 // Array of coefficients for the user-defined objective.
695 // Indexed by column number. Used in Phase-II.
696 DenseRow objective_;
697
698 // Objective offset and scaling factor of the linear program given to Solve().
699 // This is used to display the correct objective values in the logs with
700 // ComputeInitialProblemObjectiveValue().
701 Fractional objective_offset_;
702 Fractional objective_scaling_factor_;
703
704 // Used in dual phase I to keep track of the non-basic dual infeasible
705 // columns and their sign of infeasibility (+1 or -1).
706 DenseRow dual_infeasibility_improvement_direction_;
707 int num_dual_infeasible_positions_;
708
709 // A temporary scattered column that is always reset to all zero after use.
710 ScatteredColumn initially_all_zero_scratchpad_;
711
712 // Array of column index, giving the column number corresponding
713 // to a given basis row.
714 RowToColMapping basis_;
715 RowToColMapping tmp_basis_;
716
717 // Vector of strings containing the names of variables.
718 // Indexed by column number.
719 StrictITIVector<ColIndex, std::string> variable_name_;
720
721 // Only used for logging. What triggered a refactorization.
722 RefactorizationReason last_refactorization_reason_;
723
724 // Information about the solution computed by the last Solve().
725 Fractional solution_objective_value_;
726 DenseColumn solution_dual_values_;
727 DenseRow solution_reduced_costs_;
728 DenseRow solution_primal_ray_;
729 DenseColumn solution_dual_ray_;
730 DenseRow solution_dual_ray_row_combination_;
731 BasisState solution_state_;
732 bool solution_state_has_been_set_externally_;
733
734 // If this is cleared, we assume they are none.
735 DenseRow variable_starting_values_;
736
737 // See MutableTransposedMatrixWithSlack().
738 bool transpose_was_changed_ = false;
739
740 // This is known as 'd' in the literature and is set during each pivot to the
741 // right inverse of the basic entering column of A by ComputeDirection().
742 // ComputeDirection() also fills direction_.non_zeros with the position of the
743 // non-zero.
744 ScatteredColumn direction_;
745 Fractional direction_infinity_norm_;
746
747 // Used to compute the error 'b - A.x' or 'a - B.d'.
748 DenseColumn error_;
749
750 // A random number generator. In test we use absl_random_ to have a
751 // non-deterministic behavior and avoid client depending on a golden optimal
752 // solution which prevent us from easily changing the solver.
753 random_engine_t deterministic_random_;
754 absl::BitGen absl_random_;
755
756 // A reference to one of the above random generators. Fixed at construction.
757 absl::BitGenRef random_;
758
759 // Helpers for logging the solve progress.
760 SolverLogger default_logger_;
761 SolverLogger* logger_ = &default_logger_;
762
763 // Representation of matrix B using eta matrices and LU decomposition.
764 BasisFactorization basis_factorization_;
765
766 // Classes responsible for maintaining the data of the corresponding names.
767 VariablesInfo variables_info_;
768 PrimalEdgeNorms primal_edge_norms_;
769 DualEdgeNorms dual_edge_norms_;
770 DynamicMaximum<RowIndex> dual_prices_;
771 VariableValues variable_values_;
772 UpdateRow update_row_;
773 ReducedCosts reduced_costs_;
774 EnteringVariable entering_variable_;
775 PrimalPrices primal_prices_;
776
777 // Used in dual phase I to hold the price of each possible leaving choices.
778 DenseColumn dual_pricing_vector_;
779 DenseColumn tmp_dual_pricing_vector_;
780
781 // Temporary memory used by DualMinimize().
782 std::vector<ColIndex> bound_flip_candidates_;
783
784 // Total number of iterations performed.
785 uint64_t num_iterations_ = 0;
786
787 // Number of iterations performed during the first (feasibility) phase.
788 uint64_t num_feasibility_iterations_ = 0;
789
790 // Number of iterations performed during the second (optimization) phase.
791 uint64_t num_optimization_iterations_ = 0;
792
793 // Number of iterations performed during the push/crossover phase.
794 uint64_t num_push_iterations_ = 0;
795
796 // Deterministic time for DualPhaseIUpdatePriceOnReducedCostChange().
797 int64_t num_update_price_operations_ = 0;
798
799 // Total time spent in Solve().
800 double total_time_ = 0.0;
801
802 // Time spent in the first (feasibility) phase.
803 double feasibility_time_ = 0.0;
804
805 // Time spent in the second (optimization) phase.
806 double optimization_time_ = 0.0;
807
808 // Time spent in the push/crossover phase.
809 double push_time_ = 0.0;
810
811 // The internal deterministic time during the most recent call to
812 // RevisedSimplex::AdvanceDeterministicTime.
813 double last_deterministic_time_update_ = 0.0;
814
815 // Statistics about the iterations done by PrimalMinimize().
816 IterationStats iteration_stats_;
817
818 mutable RatioTestStats ratio_test_stats_;
819
820 // Placeholder for all the function timing stats.
821 // Mutable because we time const functions like ChooseLeavingVariableRow().
822 mutable StatsGroup function_stats_;
823
824 // Proto holding all the parameters of this algorithm.
825 //
826 // Note that parameters_ may actually change during a solve as the solver may
827 // dynamically adapt some values. It is why we store the argument of the last
828 // SetParameters() call in initial_parameters_ so the next Solve() can reset
829 // it correctly.
830 GlopParameters parameters_;
831 GlopParameters initial_parameters_;
832
833 // LuFactorization used to test if a pivot will cause the new basis to
834 // not be factorizable.
835 LuFactorization test_lu_;
836
837 // Number of degenerate iterations made just before the current iteration.
838 int num_consecutive_degenerate_iterations_;
839
840 // Indicate the current phase of the solve.
841 Phase phase_ = Phase::FEASIBILITY;
842
843 // Indicates whether simplex ended due to the objective limit being reached.
844 // Note that it's not enough to compare the final objective value with the
845 // limit due to numerical issues (i.e., the limit which is reached within
846 // given tolerance on the internal objective may no longer be reached when the
847 // objective scaling and offset are taken into account).
848 bool objective_limit_reached_;
849
850 // Temporary SparseColumn used by ChooseLeavingVariableRow().
851 SparseColumn leaving_candidates_;
852
853 // Temporary vector used to hold the best leaving column candidates that are
854 // tied using the current choosing criteria. We actually only store the tied
855 // candidate #2, #3, ...; because the first tied candidate is remembered
856 // anyway.
857 std::vector<RowIndex> equivalent_leaving_choices_;
858
859 // This is used by Polish().
860 DenseRow integrality_scale_;
861};
862
863// Hides the details of the dictionary matrix implementation. In the future,
864// GLOP will support generating the dictionary one row at a time without having
865// to store the whole matrix in memory.
867 public:
870 // RevisedSimplex cannot be passed const because we have to call a non-const
871 // method ComputeDictionary.
872 // TODO(user): Overload this to take RevisedSimplex* alone when the
873 // caller would normally pass a nullptr for col_scales so this and
874 // ComputeDictionary can take a const& argument.
875 RevisedSimplexDictionary(const DenseRow* col_scales,
876 RevisedSimplex* revised_simplex)
877 : dictionary_(
878 ABSL_DIE_IF_NULL(revised_simplex)->ComputeDictionary(col_scales)),
879 basis_vars_(ABSL_DIE_IF_NULL(revised_simplex)->GetBasisVector()) {}
880
881 // This type is neither copyable nor movable.
885 ConstIterator begin() const { return dictionary_.begin(); }
886 ConstIterator end() const { return dictionary_.end(); }
888 size_t NumRows() const { return dictionary_.size(); }
890 // TODO(user): This function is a better fit for the future custom iterator.
891 ColIndex GetBasicColumnForRow(RowIndex r) const { return basis_vars_[r]; }
892 SparseRow GetRow(RowIndex r) const { return dictionary_[r]; }
894 private:
895 const RowMajorSparseMatrix dictionary_;
896 const RowToColMapping basis_vars_;
897};
898
899// TODO(user): When a row-by-row generation of the dictionary is supported,
900// implement DictionaryIterator class that would call it inside operator*().
901
902} // namespace glop
903} // namespace operations_research
904
905#endif // OR_TOOLS_GLOP_REVISED_SIMPLEX_H_
Statistic on the distribution of a sequence of integers.
Definition stats.h:288
Base class to print a nice summary of a group of statistics.
Definition stats.h:128
StatsGroup(absl::string_view name)
Definition stats.h:135
RevisedSimplexDictionary(const DenseRow *col_scales, RevisedSimplex *revised_simplex)
RevisedSimplexDictionary & operator=(const RevisedSimplexDictionary &)=delete
RowMajorSparseMatrix::const_iterator ConstIterator
Entry point of the revised simplex algorithm implementation.
Fractional GetVariableValue(ColIndex col) const
void SetIntegralityScale(ColIndex col, Fractional scale)
ConstraintStatus GetConstraintStatus(RowIndex row) const
void LoadStateForNextSolve(const BasisState &state)
Uses the given state as a warm-start for the next Solve() call.
const GlopParameters & GetParameters() const
Fractional GetDualValue(RowIndex row) const
Fractional GetReducedCost(ColIndex col) const
const CompactSparseMatrix & MatrixWithSlack() const
const ScatteredRow & GetUnitRowLeftInverse(RowIndex row)
Fractional GetConstraintActivity(RowIndex row) const
const DenseBitRow & GetNotBasicBitRow() const
RowMajorSparseMatrix ComputeDictionary(const DenseRow *column_scales)
const DenseRow & GetDualRayRowCombination() const
This is the "dual ray" linear combination of the matrix rows.
RowIndex GetProblemNumRows() const
Getters to retrieve all the information computed by the last Solve().
void ComputeBasicVariablesForState(const LinearProgram &linear_program, const BasisState &state)
std::string StatString()
Returns statistics about this class as a string.
const BasisFactorization & GetBasisFactorization() const
CompactSparseMatrix * MutableTransposedMatrixWithSlack()
VariableStatus GetVariableStatus(ColIndex col) const
void SetParameters(const GlopParameters &parameters)
Sets or gets the algorithm parameters to be used on the next Solve().
ABSL_MUST_USE_RESULT Status Solve(const LinearProgram &lp, TimeLimit *time_limit)
ABSL_MUST_USE_RESULT Status MinimizeFromTransposedMatrixWithSlack(const DenseRow &objective, Fractional objective_scaling_factor, Fractional objective_offset, TimeLimit *time_limit)
void SetStartingVariableValuesForNextSolve(const DenseRow &values)
RevisedSimplex & operator=(const RevisedSimplex &)=delete
StrictITISpan< RowIndex, const Fractional > ConstView
Definition lp_types.h:291
time_limit
Definition solve.cc:22
StrictITIVector< RowIndex, ColIndex > RowToColMapping
Definition lp_types.h:394
Bitset64< ColIndex > DenseBitRow
Row of bits.
Definition lp_types.h:375
util_intops::StrongVector< RowIndex, SparseRow > RowMajorSparseMatrix
A matrix stored by rows.
Definition sparse_row.h:61
VariableType
Different types of variables.
Definition lp_types.h:178
StrictITIVector< RowIndex, Fractional > DenseColumn
Column-vector types. Column-vector types are indexed by a row index.
Definition lp_types.h:380
StrictITIVector< ColIndex, Fractional > DenseRow
Row-vector types. Row-vector types are indexed by a column index.
Definition lp_types.h:351
ProblemStatus
Different statuses for a given problem.
Definition lp_types.h:105
VectorXd ReducedCosts(const PrimalDualHybridGradientParams &params, const ShardedQuadraticProgram &sharded_qp, const VectorXd &primal_solution, const VectorXd &dual_solution, bool use_zero_primal_objective)
In SWIG mode, we don't want anything besides these top-level includes.
std::mt19937_64 random_engine_t
util_intops::StrongVector< ColumnEntryIndex, ElementIndex > SparseColumn