Coverage for picos/solvers/solver

1# ------------------------------------------------------------------------------

4# This file is part of PICOS.

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)

100

101 self._numVars = 0

102 """Total number of scalar variables passed to OSQP."""

103

104 self._osqpVarOffset = {}

105 """Maps a PICOS variable to its column in the constraint matrix."""

106

107 self._osqpConOffset = {}

108 """Maps a PICOS constraint to its row in the constraint matrix."""

109

110 self._objectiveOffset = 0.0

111 """Objective function constant offset."""

112

113 def reset_problem(self):

114 """Implement :meth:`~.solver.Solver.reset_problem`."""

115 self.int = None

116

117 self._numVars = 0

118 self._osqpVarOffset.clear()

119 self._osqpConOffset.clear()

120 self._objectiveOffset = 0.0

121

122 def _affine_expression_to_G_and_h(self, expression):

123 assert isinstance(expression, AffineExpression)

124

125 return expression.scipy_sparse_matrix_form(

126 varOffsetMap=self._osqpVarOffset, dense_b=True)

127

128 _Gh = _affine_expression_to_G_and_h

129

130 def _import_variables(self):

131 offset = 0

132 for variable in self.ext.variables.values():

133 dim = variable.dim

134

135 # Register the variable.

136 self._osqpVarOffset[variable] = offset

137 offset += dim

138

139 assert offset == self._numVars

140

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)

147

148 def _import_objective(self):

149 direction, objective = self.ext.no

150

151 # OSQP only supports minimization; flip the sign for maximization.

152 if direction == "max":

153 objective = -objective

154

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."

164

165 # Import quadratic part.

166 for xy, Q in sparse_quads.items():

167 x, y = xy

168 m, n = x.dim, y.dim

169

170 dx = self._osqpVarOffset[x]

171 dy = self._osqpVarOffset[y]

172

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

178

179 # OSQP adds a factor of 0.5; cancel it.

180 Q = 2*Q

181

182 # Convert from cvxopt sparse to scipy sparse matrix.

183 Q = cvx2csc(Q)

184

185 self.int["P"][dx:dx + m, dy:dy + n] = Q

186

187 # Import linear part.

188 q, c = self._Gh(affine_part)

189

190 self.int["q"] = numpy.ravel(q.todense())

191 self._objectiveOffset = float(c[0])

192

193 def _import_affine_constraint(self, constraint):

194 from scipy import sparse

195

196 assert isinstance(constraint, AffineConstraint)

197

198 A, minus_b = self._Gh(constraint.lmr)

199 b = -minus_b

200

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"))

209

210 self._osqpConOffset[constraint] = len(self.int["l"])

211

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])

216

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__)

224

225 def _import_problem(self):

226 from scipy import sparse

227

228 self._numVars = n = sum(var.dim for var in self.ext.variables.values())

229

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)),

239

240 # Objective function linear coefficients.

241 "q": numpy.zeros(n),

242

243 # Linear inequality coefficient matrix.

244 "A": sparse.csc_matrix((0, n)),

245

246 # Linear inequality lower bound.

247 "l": numpy.zeros(0),

248

249 # Linear inequality upper bound.

250 "u": numpy.zeros(0),

251 }

252

253 # Import variables without their bounds.

254 self._import_variables()

255

256 # Set objective.

257 self._import_objective()

258

259 # Import constraints.

260 for constraint in self.ext.constraints.values():

261 self._import_constraint(constraint)

262

263 # Convert from LIL to CSC manually to avoid a warning by OSQP.

264 self.int["P"] = self.int["P"].tocsc()

265

266 def _update_problem(self):

267 raise NotImplementedError

268

269 def _solve(self):

270 import osqp

271

272 options = {}

273

274 # verbosity

275 options["verbose"] = (self.verbosity() >= 1)

276

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

280

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

284

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

288

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

292

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)

301

302 # timelimit

303 if self.ext.options.timelimit is not None:

304 options["time_limit"] = self.ext.options.timelimit

305

306 # Enable polishing to increase chance of obeying precision limits.

307 options["polish"] = True

308

309 # Handle OSQP-specific options.

310 options.update(self.ext.options.osqp_params)

311

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")

316

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)

322

323 model = osqp.OSQP()

324 model.setup(**self.int, **options)

325

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

333

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])

340

341 if None in primal:

342 primal = None

343 else:

344 primal = cvxopt.matrix(primal)

345

346 primals[variable] = primal

347

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

356

357 assert isinstance(constraint, AffineConstraint)

358

359 offset = self._osqpConOffset[constraint]

360 length = len(constraint)

361 dual = list(result.y[offset:offset + length])

362

363 if None in dual:

364 dual = None

365 else:

366 dual = cvxopt.matrix(dual)

367

368 if not constraint.is_increasing():

369 dual = -dual

370

371 duals[constraint] = dual

372

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

378

379 # Flip back the sign for maximization.

380 if self.ext.no.direction == "max":

381 value = -value

382

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)

407

408 return self._make_solution(

409 value, primals, duals, primalStatus, dualStatus, problemStatus,

410 {"osqp_info": result.info if result else None})

411

412

413# --------------------------------------

414__all__ = api_end(_API_START, globals())

Coverage for picos/solvers/solver_osqp.py: 89.57%

211 statements