Coverage for lpsolvers / pdlp_.py: 82%

1#!/usr/bin/env python 

2# -*- coding: utf-8 -*- 

3# 

4# SPDX-License-Identifier: LGPL-3.0-or-later 

5# Copyright 2023 Inria 

6 

7"""Solver interface for `PDLP`_. 

8 

9.. _PDLP: https://developers.google.com/optimization/lp/pdlp_math 

10 

11PDLP is a first-order method for convex quadratic programming aiming for 

12high-accuracy solutions and scaling to large problems. If you use PDLP in your 

13academic works, consider citing the corresponding paper [Applegate2021]_. 

14""" 

15 

16from typing import Optional 

17 

18import numpy as np 

19import scipy.sparse as spa 

20from ortools.pdlp import solve_log_pb2, solvers_pb2 

21from ortools.pdlp.python import pdlp 

22 

23 

24def pdlp_solve_lp( 

25 c: np.ndarray, 

26 G: np.ndarray, 

27 h: np.ndarray, 

28 A: Optional[np.ndarray] = None, 

29 b: Optional[np.ndarray] = None, 

30 verbose: bool = False, 

31 eps_optimal_absolute: Optional[float] = None, 

32 eps_optimal_relative: Optional[float] = None, 

33 time_sec_limits: Optional[float] = None, 

34 **kwargs, 

35) -> np.ndarray: 

36 """Solve a quadratic program using PDLP. 

37 

38 Parameters 

39 ---------- 

40 c : 

41 Linear cost vector. 

42 G : 

43 Linear inequality constraint matrix. 

44 h : 

45 Linear inequality constraint vector. 

46 A : 

47 Linear equality constraint matrix. 

48 b : 

49 Linear equality constraint vector. 

50 verbose : 

51 Set to `True` to print out extra information. 

52 verbose : 

53 Set to `True` to print out extra information. 

54 eps_optimal_absolute : 

55 Absolute tolerance on the primal-dual residuals and duality gap. See 

56 *e.g.* [tolerances]_ for an overview of solver tolerances. 

57 eps_optimal_relative : 

58 Relative tolerance on the primal-dual residuals and duality gap. See 

59 *e.g.* [tolerances]_ for an overview of solver tolerances. 

60 time_sec_limits : 

61 Maximum computation time the solver is allowed, in seconds. 

62 

63 Returns 

64 ------- 

65 : 

66 Primal solution to the QP, if found, otherwise ``None``. 

67 

68 Notes 

69 ----- 

70 All other keyword arguments are forwarded as parameters to PDLP. For 

71 instance, you can call ``pdlp_solve_qp(P, q, G, h, num_threads=3, 

72 verbosity_level=2)``. For a quick overview, the solver accepts the 

73 following settings: 

74 

75 .. list-table:: 

76 :widths: 30 70 

77 :header-rows: 1 

78 

79 * - Name 

80 - Effect 

81 * - ``num_threads`` 

82 - Number of threads to use (positive). 

83 * - ``verbosity_level`` 

84 - Verbosity level from 0 (no logging) to 4 (extensive logging). 

85 * - ``initial_primal_weight`` 

86 - Initial value of the primal weight (ratio of primal over dual step 

87 sizes). 

88 * - ``l_inf_ruiz_iterations`` 

89 - Number of L-infinity Ruiz rescaling iterations applied to the 

90 constraint matrix. 

91 * - ``l2_norm_rescaling`` 

92 - If set to ``True``, applies L2-norm rescaling after Ruiz rescaling. 

93 

94 This list is not exhaustive. Check out the solver's `Protocol Bufffers file 

95 <https://github.com/google/or-tools/blob/8768ed7a43f8899848effb71295a790f3ecbe2f2/ortools/pdlp/solvers.proto>`__ 

96 for more. See also the `Mathematical background for PDLP 

97 <https://developers.google.com/optimization/lp/pdlp_math>`__. 

98 """ 

99 n = c.shape[0] 

100 

101 A_pdlp = None 

102 lc_pdlp = None 

103 uc_pdlp = None 

104 if G is not None and h is not None: 

105 A_pdlp = G 

106 lc_pdlp = np.full(h.shape, -np.inf) 

107 uc_pdlp = h 

108 if A is not None and b is not None: 

109 A_pdlp = A if A_pdlp is None else spa.vstack([A_pdlp, A], format="csc") 

110 lc_pdlp = b if lc_pdlp is None else np.hstack([lc_pdlp, b]) 

111 uc_pdlp = b if uc_pdlp is None else np.hstack([uc_pdlp, b]) 

112 lv_pdlp = np.full((n,), -np.inf) # custom lb vector can go here 

113 uv_pdlp = np.full((n,), +np.inf) # custom ub vector can go here 

114 

115 qp = pdlp.QuadraticProgram() 

116 # qp.objective_matrix = np.diag(...) 

117 qp.objective_vector = c 

118 if A_pdlp is not None: 

119 qp.constraint_matrix = A_pdlp 

120 qp.constraint_lower_bounds = lc_pdlp 

121 qp.constraint_upper_bounds = uc_pdlp 

122 qp.variable_lower_bounds = lv_pdlp 

123 qp.variable_upper_bounds = uv_pdlp 

124 

125 params = solvers_pb2.PrimalDualHybridGradientParams() 

126 optimality = params.termination_criteria.simple_optimality_criteria 

127 if eps_optimal_absolute is not None: 

128 optimality.eps_optimal_absolute = eps_optimal_absolute 

129 if eps_optimal_relative is not None: 

130 optimality.eps_optimal_relative = eps_optimal_relative 

131 if time_sec_limits is not None: 

132 params.termination_criteria.time_sec_limits = time_sec_limits 

133 if verbose and "verbosity_level" not in kwargs: 

134 params.verbosity_level = 1 if verbose else 0 

135 for param, value in kwargs.items(): 

136 setattr(params, param, value) 

137 

138 result = pdlp.primal_dual_hybrid_gradient(qp, params) 

139 log = result.solve_log 

140 if log.termination_reason != solve_log_pb2.TERMINATION_REASON_OPTIMAL: 

141 raise ValueError("Linear program is not feasible") 

142 return result.primal_solution