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

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# ------------------------------------------------------------------------------ 

18 

19"""Implementation of :class:`OSQPSolver`.""" 

20 

21import cvxopt 

22import numpy 

23 

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 

34 

35_API_START = api_start(globals()) 

36# ------------------------------- 

37 

38 

39class OSQPSolver(Solver): 

40 """Interface to the OSQP solver.""" 

41 

42 SUPPORTED = Specification( 

43 objectives=[ 

44 AffineExpression, 

45 QuadraticExpression], 

46 variables=CONTINUOUS_VARTYPES, 

47 constraints=[ 

48 DummyConstraint, 

49 AffineConstraint]) 

50 

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 

57 

58 if footprint.nonconvex_quadratic_objective: 

59 if explain: 

60 return (False, "QPs with a nonconvex objective.") 

61 else: 

62 return False 

63 

64 if footprint not in cls.SUPPORTED: 

65 if explain: 

66 return False, cls.SUPPORTED.mismatch_reason(footprint) 

67 else: 

68 return False 

69 

70 return (True, None) if explain else True 

71 

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 

78 

79 @classmethod 

80 def test_availability(cls): 

81 """Implement :meth:`~.solver.Solver.test_availability`.""" 

82 cls.check_import("osqp") 

83 

84 @classmethod 

85 def names(cls): 

86 """Implement :meth:`~.solver.Solver.names`.""" 

87 return "osqp", "OSQP", "Operator Splitting QP Solver", None 

88 

89 @classmethod 

90 def is_free(cls): 

91 """Implement :meth:`~.solver.Solver.is_free`.""" 

92 return True 

93 

94 def __init__(self, problem): 

95 """Initialize a OSQP solver interface. 

96 

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