Coverage for picos/solvers/solver_osqp.py: 89.57%
211 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
1# ------------------------------------------------------------------------------
2# Copyright (C) 2021-2022 Maximilian Stahlberg
3#
4# This file is part of PICOS.
5#
6# PICOS is free software: you can redistribute it and/or modify it under the
7# terms of the GNU General Public License as published by the Free Software
8# Foundation, either version 3 of the License, or (at your option) any later
9# version.
10#
11# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
12# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License along with
16# this program. If not, see <http://www.gnu.org/licenses/>.
17# ------------------------------------------------------------------------------
19"""Implementation of :class:`OSQPSolver`."""
21import cvxopt
22import numpy
24from ..apidoc import api_end, api_start
25from ..constraints import AffineConstraint, DummyConstraint
26from ..expressions import (CONTINUOUS_VARTYPES, AffineExpression,
27 QuadraticExpression)
28from ..expressions.data import cvx2csc
29from ..modeling.footprint import Specification
30from ..modeling.solution import (PS_FEASIBLE, PS_INFEASIBLE, PS_UNBOUNDED,
31 PS_UNKNOWN, SS_FAILURE, SS_INFEASIBLE,
32 SS_OPTIMAL, SS_PREMATURE, SS_UNKNOWN)
33from .solver import Solver
35_API_START = api_start(globals())
36# -------------------------------
39class OSQPSolver(Solver):
40 """Interface to the OSQP solver."""
42 SUPPORTED = Specification(
43 objectives=[
44 AffineExpression,
45 QuadraticExpression],
46 variables=CONTINUOUS_VARTYPES,
47 constraints=[
48 DummyConstraint,
49 AffineConstraint])
51 @classmethod
52 def supports(cls, footprint, explain=False):
53 """Implement :meth:`~.solver.Solver.supports`."""
54 result = Solver.supports(footprint, explain)
55 if not result or (explain and not result[0]):
56 return result
58 if footprint.nonconvex_quadratic_objective:
59 if explain:
60 return (False, "QPs with a nonconvex objective.")
61 else:
62 return False
64 if footprint not in cls.SUPPORTED:
65 if explain:
66 return False, cls.SUPPORTED.mismatch_reason(footprint)
67 else:
68 return False
70 return (True, None) if explain else True
72 @classmethod
73 def default_penalty(cls):
74 """Implement :meth:`~.solver.Solver.default_penalty`."""
75 # OSQP is an established free/open source solver but it has issues with
76 # moderate numbers of affine inequalities, so leave it to user choice.
77 return 2.0
79 @classmethod
80 def test_availability(cls):
81 """Implement :meth:`~.solver.Solver.test_availability`."""
82 cls.check_import("osqp")
84 @classmethod
85 def names(cls):
86 """Implement :meth:`~.solver.Solver.names`."""
87 return "osqp", "OSQP", "Operator Splitting QP Solver", None
89 @classmethod
90 def is_free(cls):
91 """Implement :meth:`~.solver.Solver.is_free`."""
92 return True
94 def __init__(self, problem):
95 """Initialize a OSQP solver interface.
97 :param ~picos.Problem problem: The problem to be solved.
98 """
99 super(OSQPSolver, self).__init__(problem)
101 self._numVars = 0
102 """Total number of scalar variables passed to OSQP."""
104 self._osqpVarOffset = {}
105 """Maps a PICOS variable to its column in the constraint matrix."""
107 self._osqpConOffset = {}
108 """Maps a PICOS constraint to its row in the constraint matrix."""
110 self._objectiveOffset = 0.0
111 """Objective function constant offset."""
113 def reset_problem(self):
114 """Implement :meth:`~.solver.Solver.reset_problem`."""
115 self.int = None
117 self._numVars = 0
118 self._osqpVarOffset.clear()
119 self._osqpConOffset.clear()
120 self._objectiveOffset = 0.0
122 def _affine_expression_to_G_and_h(self, expression):
123 assert isinstance(expression, AffineExpression)
125 return expression.scipy_sparse_matrix_form(
126 varOffsetMap=self._osqpVarOffset, dense_b=True)
128 _Gh = _affine_expression_to_G_and_h
130 def _import_variables(self):
131 offset = 0
132 for variable in self.ext.variables.values():
133 dim = variable.dim
135 # Register the variable.
136 self._osqpVarOffset[variable] = offset
137 offset += dim
139 assert offset == self._numVars
141 # Add variable bounds as affine constraints.
142 for variable in self.ext.variables.values():
143 # TODO: Import lower and upper bound in a single constraint instead.
144 bounds = variable.bound_constraint
145 if bounds:
146 self._import_affine_constraint(bounds)
148 def _import_objective(self):
149 direction, objective = self.ext.no
151 # OSQP only supports minimization; flip the sign for maximization.
152 if direction == "max":
153 objective = -objective
155 # Split objective into quadratic and affine part.
156 if isinstance(objective, AffineExpression):
157 sparse_quads = {}
158 affine_part = objective
159 elif isinstance(objective, QuadraticExpression):
160 sparse_quads = objective._sparse_quads
161 affine_part = objective.aff
162 else:
163 assert False, "Unexpected objective."
165 # Import quadratic part.
166 for xy, Q in sparse_quads.items():
167 x, y = xy
168 m, n = x.dim, y.dim
170 dx = self._osqpVarOffset[x]
171 dy = self._osqpVarOffset[y]
173 # OSQP reads only the upper triangular part.
174 if dx > dy:
175 dx, dy = dy, dx
176 m, n = n, m
177 Q = Q.T
179 # OSQP adds a factor of 0.5; cancel it.
180 Q = 2*Q
182 # Convert from cvxopt sparse to scipy sparse matrix.
183 Q = cvx2csc(Q)
185 self.int["P"][dx:dx + m, dy:dy + n] = Q
187 # Import linear part.
188 q, c = self._Gh(affine_part)
190 self.int["q"] = numpy.ravel(q.todense())
191 self._objectiveOffset = float(c[0])
193 def _import_affine_constraint(self, constraint):
194 from scipy import sparse
196 assert isinstance(constraint, AffineConstraint)
198 A, minus_b = self._Gh(constraint.lmr)
199 b = -minus_b
201 if constraint.is_equality():
202 l = u = b
203 elif constraint.is_increasing():
204 l = numpy.full(len(b), float("-inf"))
205 u = b
206 elif constraint.is_decreasing():
207 l = b
208 u = numpy.full(len(b), float("+inf"))
210 self._osqpConOffset[constraint] = len(self.int["l"])
212 # TODO: Add matrices to a list and concatenate them at once later.
213 self.int["A"] = sparse.vstack([self.int["A"], A], format="csc")
214 self.int["l"] = numpy.concatenate([self.int["l"], l])
215 self.int["u"] = numpy.concatenate([self.int["u"], u])
217 def _import_constraint(self, constraint):
218 if isinstance(constraint, AffineConstraint):
219 self._import_affine_constraint(constraint)
220 else:
221 assert isinstance(constraint, DummyConstraint), \
222 "Unexpected constraint type: {}".format(
223 constraint.__class__.__name__)
225 def _import_problem(self):
226 from scipy import sparse
228 self._numVars = n = sum(var.dim for var in self.ext.variables.values())
230 # OSQP's internal problem representation is stateful but supports
231 # updates only by setting or replacing whole matrices and not on a
232 # per-variable or per-constraint basis. We thus pretend it was stateless
233 # and use the osqp.solve function in a similar manner as with CVXOPT.
234 # TODO: Consider supporting the limited update capabilities of OSQP.
235 self.int = {
236 # Objective function quadratic form.
237 # NOTE: lil_matrix for cheap updates, converted to csc_matrix later.
238 "P": sparse.lil_matrix((n, n)),
240 # Objective function linear coefficients.
241 "q": numpy.zeros(n),
243 # Linear inequality coefficient matrix.
244 "A": sparse.csc_matrix((0, n)),
246 # Linear inequality lower bound.
247 "l": numpy.zeros(0),
249 # Linear inequality upper bound.
250 "u": numpy.zeros(0),
251 }
253 # Import variables without their bounds.
254 self._import_variables()
256 # Set objective.
257 self._import_objective()
259 # Import constraints.
260 for constraint in self.ext.constraints.values():
261 self._import_constraint(constraint)
263 # Convert from LIL to CSC manually to avoid a warning by OSQP.
264 self.int["P"] = self.int["P"].tocsc()
266 def _update_problem(self):
267 raise NotImplementedError
269 def _solve(self):
270 import osqp
272 options = {}
274 # verbosity
275 options["verbose"] = (self.verbosity() >= 1)
277 # abs_prim_fsb_tol
278 if self.ext.options.abs_prim_fsb_tol is not None:
279 options["eps_prim_inf"] = self.ext.options.abs_prim_fsb_tol
281 # abs_dual_fsb_tol
282 if self.ext.options.abs_dual_fsb_tol is not None:
283 options["eps_dual_inf"] = self.ext.options.abs_dual_fsb_tol
285 # abs_ipm_opt_tol
286 if self.ext.options.abs_ipm_opt_tol is not None:
287 options["eps_abs"] = self.ext.options.abs_ipm_opt_tol
289 # rel_ipm_opt_tol
290 if self.ext.options.rel_ipm_opt_tol is not None:
291 options["eps_rel"] = self.ext.options.rel_ipm_opt_tol
293 # max_iterations
294 # NOTE: OSQP can hang long already for moderately sized LPs but we still
295 # "remove" the iteration limit to obey the PICOS setting. This is
296 # fine as long as OSQP requires user selection.
297 if self.ext.options.max_iterations is not None:
298 options["max_iter"] = self.ext.options.max_iterations
299 else:
300 options["max_iter"] = int(1e9)
302 # timelimit
303 if self.ext.options.timelimit is not None:
304 options["time_limit"] = self.ext.options.timelimit
306 # Enable polishing to increase chance of obeying precision limits.
307 options["polish"] = True
309 # Handle OSQP-specific options.
310 options.update(self.ext.options.osqp_params)
312 # Handle unsupported options.
313 # TODO: Support hotstart.
314 self._handle_unsupported_options("lp_root_method", "lp_node_method",
315 "treememory", "max_fsb_nodes", "hotstart")
317 # Attempt to solve the problem.
318 with self._header(), self._stopwatch():
319 # NOTE: There is supposed to be a direct function osqp.solve but it
320 # does not exist for my installation of 0.6.2.
321 # result = osqp.solve(**self.int, **options)
323 model = osqp.OSQP()
324 model.setup(**self.int, **options)
326 try:
327 result = model.solve()
328 except ValueError as error:
329 if str(error) == "OSQP solve error!":
330 result = None
331 else:
332 raise
334 # Retrieve primals.
335 primals = {}
336 if result and self.ext.options.primals is not False:
337 for variable in self.ext.variables.values():
338 offset = self._osqpVarOffset[variable]
339 primal = list(result.x[offset:offset + variable.dim])
341 if None in primal:
342 primal = None
343 else:
344 primal = cvxopt.matrix(primal)
346 primals[variable] = primal
348 # Retrieve duals.
349 duals = {}
350 if result and self.ext.options.duals is not False:
351 for constraint in self.ext.constraints.values():
352 if isinstance(constraint, DummyConstraint):
353 duals[constraint] = cvxopt.spmatrix(
354 [], [], [], constraint.size)
355 continue
357 assert isinstance(constraint, AffineConstraint)
359 offset = self._osqpConOffset[constraint]
360 length = len(constraint)
361 dual = list(result.y[offset:offset + length])
363 if None in dual:
364 dual = None
365 else:
366 dual = cvxopt.matrix(dual)
368 if not constraint.is_increasing():
369 dual = -dual
371 duals[constraint] = dual
373 # Retrieve objective value.
374 value = result.info.obj_val if result else None
375 if value is not None:
376 # Add back the constant part.
377 value += self._objectiveOffset
379 # Flip back the sign for maximization.
380 if self.ext.no.direction == "max":
381 value = -value
383 # Retrieve solution status.
384 status = result.info.status if result else None
385 if status is None:
386 primalStatus = SS_FAILURE
387 dualStatus = SS_FAILURE
388 problemStatus = PS_UNKNOWN
389 elif status == "solved":
390 primalStatus = SS_OPTIMAL
391 dualStatus = SS_OPTIMAL
392 problemStatus = PS_FEASIBLE
393 elif status == "primal infeasible":
394 primalStatus = SS_INFEASIBLE
395 dualStatus = SS_UNKNOWN
396 problemStatus = PS_INFEASIBLE
397 elif status == "dual infeasible":
398 primalStatus = SS_UNKNOWN
399 dualStatus = SS_INFEASIBLE
400 problemStatus = PS_UNBOUNDED
401 elif status in ("maximum iterations reached", "run time limit reached"):
402 primalStatus = SS_PREMATURE
403 dualStatus = SS_PREMATURE
404 problemStatus = PS_UNKNOWN
405 else:
406 assert False, "Unknown solver status '{}'".format(status)
408 return self._make_solution(
409 value, primals, duals, primalStatus, dualStatus, problemStatus,
410 {"osqp_info": result.info if result else None})
413# --------------------------------------
414__all__ = api_end(_API_START, globals())