Coverage for picos/solvers/solver_ecos.py: 83.78%

296 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-03-26 07:46 +0000

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

2# Copyright (C) 2018-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:`ECOSSolver`.""" 

20 

21import cvxopt 

22import numpy 

23 

24from ..apidoc import api_end, api_start 

25from ..constraints import (AffineConstraint, DummyConstraint, 

26 ExpConeConstraint, RSOCConstraint, SOCConstraint) 

27from ..expressions import AffineExpression, BinaryVariable, IntegerVariable 

28from ..modeling.footprint import Specification 

29from ..modeling.solution import (PS_FEASIBLE, PS_INFEASIBLE, PS_UNBOUNDED, 

30 PS_UNKNOWN, PS_UNSTABLE, SS_FAILURE, 

31 SS_INFEASIBLE, SS_OPTIMAL, SS_PREMATURE, 

32 SS_UNKNOWN) 

33from .solver import Solver 

34 

35_API_START = api_start(globals()) 

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

37 

38 

39class ECOSSolver(Solver): 

40 """Interface to the ECOS solver via its official Python interface.""" 

41 

42 SUPPORTED = Specification( 

43 objectives=[ 

44 AffineExpression], 

45 constraints=[ 

46 DummyConstraint, 

47 AffineConstraint, 

48 SOCConstraint, 

49 RSOCConstraint, 

50 ExpConeConstraint]) 

51 

52 @classmethod 

53 def supports(cls, footprint, explain=False): 

54 """Implement :meth:`~.solver.Solver.supports`.""" 

55 result = Solver.supports(footprint, explain) 

56 if not result or (explain and not result[0]): 

57 return result 

58 

59 if footprint not in cls.SUPPORTED: 

60 if explain: 

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

62 else: 

63 return False 

64 

65 return (True, None) if explain else True 

66 

67 @classmethod 

68 def default_penalty(cls): 

69 """Implement :meth:`~.solver.Solver.default_penalty`.""" 

70 return 1.0 # Stable free/open source solver. 

71 

72 @classmethod 

73 def test_availability(cls): 

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

75 cls.check_import("ecos") 

76 

77 @classmethod 

78 def names(cls): 

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

80 return "ecos", "ECOS", "Embedded Conic Solver", "ecos-python" 

81 

82 @classmethod 

83 def is_free(cls): 

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

85 return True 

86 

87 def __init__(self, problem): 

88 """Initialize an ECOS solver interface. 

89 

90 :param ~picos.Problem problem: The problem to be solved. 

91 """ 

92 super(ECOSSolver, self).__init__(problem) 

93 

94 self._numVars = 0 

95 """Total number of scalar variables passed to ECOS.""" 

96 

97 self._ecosVarOffset = {} 

98 """Maps a PICOS variable to its offset in the ECOS constraint matrix.""" 

99 

100 self._ecosConIndices = dict() 

101 """Maps a PICOS constraints to a range of ECOS constraint indices. 

102 

103 Note that equality constraints use a different index space.""" 

104 

105 self._ecosModuleCache = None 

106 """A cached reference to the ECOS module.""" 

107 

108 def reset_problem(self): 

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

110 self.int = None 

111 

112 self._numVars = 0 

113 self._ecosVarOffset.clear() 

114 self._ecosConIndices.clear() 

115 

116 @property 

117 def ecos(self): 

118 """Return the ECOS core module ecos.py. 

119 

120 The module is obtained by ``import ecos`` up to ECOS 2.0.6 and by 

121 ``import ecos.ecos`` starting with ECOS 2.0.7. 

122 """ 

123 if self._ecosModuleCache is None: 

124 import ecos 

125 if hasattr(ecos, "ecos"): 

126 self._ecosModuleCache = ecos.ecos 

127 else: 

128 self._ecosModuleCache = ecos 

129 return self._ecosModuleCache 

130 

131 @staticmethod 

132 def stack(*args): 

133 """Stack vectors or matrices, the latter vertically.""" 

134 import scipy.sparse 

135 

136 if isinstance(args[0], scipy.sparse.spmatrix): 

137 for i in range(1, len(args)): 

138 assert isinstance(args[i], scipy.sparse.spmatrix) 

139 

140 return scipy.sparse.vstack(args, format="csc") 

141 else: 

142 for i in range(len(args)): 

143 assert isinstance(args[i], numpy.ndarray) 

144 

145 return numpy.concatenate(args) 

146 

147 def _affine_expression_to_G_and_h(self, expression): 

148 assert isinstance(expression, AffineExpression) 

149 

150 return expression.scipy_sparse_matrix_form( 

151 varOffsetMap=self._ecosVarOffset, dense_b=True) 

152 

153 _Gh = _affine_expression_to_G_and_h 

154 

155 def _import_variables(self): 

156 offset = 0 

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

158 dim = variable.dim 

159 

160 # Mark integer variables. 

161 if isinstance(variable, (BinaryVariable, IntegerVariable)): 

162 ecosIndices = range(offset, offset + dim) 

163 

164 if isinstance(variable, BinaryVariable): 

165 self.int["bool_vars_idx"].extend(ecosIndices) 

166 else: 

167 self.int["int_vars_idx"].extend(ecosIndices) 

168 

169 # Register the variable. 

170 self._ecosVarOffset[variable] = offset 

171 offset += dim 

172 

173 # Import bounds. 

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

175 bounds = variable.bound_constraint 

176 if bounds: 

177 self._import_affine_constraint(bounds) 

178 

179 assert offset == self._numVars 

180 

181 def _append_equality(self, A, b): 

182 offset = len(self.int["b"]) 

183 length = len(b) 

184 

185 self.int["A"] = self.stack(self.int["A"], A) 

186 self.int["b"] = self.stack(self.int["b"], b) 

187 

188 return range(offset, offset + length) 

189 

190 def _append_inequality(self, G, h, typecode): 

191 assert typecode in self.int["dims"] 

192 

193 # Make sure constraints are appened in the proper order. 

194 if typecode == "l": 

195 assert len(self.int["dims"]["q"]) == 0 \ 

196 and self.int["dims"]["e"] == 0 

197 elif typecode == "q": 

198 assert self.int["dims"]["e"] == 0 

199 

200 offset = len(self.int["h"]) 

201 length = len(h) 

202 

203 self.int["G"] = self.stack(self.int["G"], G) 

204 self.int["h"] = self.stack(self.int["h"], h) 

205 

206 if typecode == "q": 

207 self.int["dims"][typecode].append(length) 

208 elif typecode == "e": 

209 self.int["dims"][typecode] += 1 

210 else: 

211 self.int["dims"][typecode] += length 

212 

213 return range(offset, offset + length) 

214 

215 def _import_affine_constraint(self, constraint): 

216 assert isinstance(constraint, AffineConstraint) 

217 

218 G, minus_h = self._Gh(constraint.le0) 

219 h = -minus_h 

220 

221 if constraint.is_equality(): 

222 return self._append_equality(G, h) 

223 else: 

224 return self._append_inequality(G, h, "l") 

225 

226 def _import_socone_constraint(self, constraint): 

227 assert isinstance(constraint, SOCConstraint) 

228 

229 (A, b) = self._Gh(constraint.ne) 

230 (c, d) = self._Gh(constraint.ub) 

231 

232 return self._append_inequality( 

233 self.stack(-c, -A), self.stack(d, b), "q") 

234 

235 def _import_rscone_constraint(self, constraint): 

236 assert isinstance(constraint, RSOCConstraint) 

237 

238 (A, b) = self._Gh(constraint.ne) 

239 (c1, d1) = self._Gh(constraint.ub1) 

240 (c2, d2) = self._Gh(constraint.ub2) 

241 

242 return self._append_inequality( 

243 self.stack(-(c1 + c2), -2.0*A, c2 - c1), 

244 self.stack(d1 + d2, 2.0*b, d1 - d2), "q") 

245 

246 def _import_expcone_constraint(self, constraint): 

247 assert isinstance(constraint, ExpConeConstraint) 

248 

249 (Gx, hx) = self._Gh(constraint.x) 

250 (Gy, hy) = self._Gh(constraint.y) 

251 (Gz, hz) = self._Gh(constraint.z) 

252 

253 # ECOS' exponential cone is cl{(x,y,z) | exp(x/z) <= y/z, z > 0}, 

254 # PICOS' is cl{(x,y,z) | x >= y*exp(z/y), y > 0}. Note that given y > 0 

255 # it is x >= y*exp(z/y) if and only if exp(z/y) <= x/y. Therefor we can 

256 # transform from our coordinates to theirs with the mapping 

257 # (x, y, z) ↦ (z, x, y). Further, G and h with G = (Gx, Gy, Gz) and 

258 # h = (hx, hy, hz) are such that G*X + h = (x, y, z) where X is the 

259 # row-vectorization of all variables. ECOS however expects G and h such 

260 # that h - G*X is constrained to be in the exponential cone. 

261 return self._append_inequality( 

262 -self.stack(Gz, Gx, Gy), self.stack(hz, hx, hy), "e") 

263 

264 def _import_objective(self): 

265 direction, objective = self.ext.no 

266 

267 # ECOS only supports minimization; flip the sign for maximization. 

268 sign = 1.0 if direction == "min" else -1.0 

269 

270 # Import coefficients. 

271 for variable, coefficient in objective._linear_coefs.items(): 

272 for localIndex in range(variable.dim): 

273 index = self._ecosVarOffset[variable] + localIndex 

274 self.int["c"][index] = sign * coefficient[localIndex] 

275 

276 def _import_problem(self): 

277 import scipy.sparse 

278 

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

280 

281 # ECOS' internal problem representation is stateless; a number of 

282 # vectors and matrices is supplied each time a search is started. 

283 # These vectors and matrices are thus stored in self.int. 

284 self.int = { 

285 # Objective function coefficients. 

286 "c": numpy.zeros(self._numVars), 

287 

288 # Linear equality left hand side. 

289 "A": scipy.sparse.csc_matrix((0, self._numVars)), 

290 

291 # Linear equality right hand side. 

292 "b": numpy.zeros(0), 

293 

294 # Conic inequality left hand side. 

295 "G": scipy.sparse.csc_matrix((0, self._numVars)), 

296 

297 # Conic inequality right hand side. 

298 "h": numpy.zeros(0), 

299 

300 # Cone definition: Linear, second order, exponential dimensions. 

301 "dims": {"l": 0, "q": [], "e": 0}, 

302 

303 # Boolean variable indices. 

304 "bool_vars_idx": [], 

305 

306 # Integer variable indices. 

307 "int_vars_idx": [] 

308 } 

309 

310 # NOTE: Constraints need to be append ordered by type. 

311 # TODO: Add fast constraints-by-type iterators to Problem. 

312 

313 # Import variables with their bounds as affine constraints. 

314 self._import_variables() 

315 

316 # Import affine constraints. 

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

318 if isinstance(constraint, AffineConstraint): 

319 self._ecosConIndices[constraint] = \ 

320 self._import_affine_constraint(constraint) 

321 

322 # Import second order cone constraints. 

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

324 if isinstance(constraint, SOCConstraint): 

325 self._ecosConIndices[constraint] = \ 

326 self._import_socone_constraint(constraint) 

327 

328 # Import rotated second order cone constraints. 

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

330 if isinstance(constraint, RSOCConstraint): 

331 self._ecosConIndices[constraint] = \ 

332 self._import_rscone_constraint(constraint) 

333 

334 # Import exponential cone constraints. 

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

336 if isinstance(constraint, ExpConeConstraint): 

337 self._ecosConIndices[constraint] = \ 

338 self._import_expcone_constraint(constraint) 

339 

340 # Make sure that no unsupported constraints are present. 

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

342 assert isinstance(constraint, (DummyConstraint, AffineConstraint, 

343 SOCConstraint, RSOCConstraint, ExpConeConstraint)) 

344 

345 # Set objective. 

346 self._import_objective() 

347 

348 def _update_problem(self): 

349 raise NotImplementedError 

350 

351 def _solve(self): 

352 ecosOptions = {} 

353 

354 # verbosity 

355 beVerbose = self.verbosity() >= 1 

356 ecosOptions["verbose"] = beVerbose 

357 ecosOptions["mi_verbose"] = beVerbose 

358 

359 # rel_prim_fsb_tol, rel_dual_fsb_tol 

360 feasibilityTols = [tol for tol in (self.ext.options.rel_prim_fsb_tol, 

361 self.ext.options.rel_dual_fsb_tol) if tol is not None] 

362 if feasibilityTols: 

363 ecosOptions["feastol"] = min(feasibilityTols) 

364 

365 # abs_ipm_opt_tol 

366 if self.ext.options.abs_ipm_opt_tol is not None: 

367 ecosOptions["abstol"] = self.ext.options.abs_ipm_opt_tol 

368 

369 # rel_ipm_opt_tol 

370 if self.ext.options.rel_ipm_opt_tol is not None: 

371 ecosOptions["reltol"] = self.ext.options.rel_ipm_opt_tol 

372 

373 # abs_bnb_opt_tol 

374 if self.ext.options.abs_bnb_opt_tol is not None: 

375 ecosOptions["mi_abs_eps"] = self.ext.options.abs_bnb_opt_tol 

376 

377 # rel_bnb_opt_tol 

378 if self.ext.options.rel_bnb_opt_tol is not None: 

379 ecosOptions["mi_rel_eps"] = self.ext.options.rel_bnb_opt_tol 

380 

381 # integrality_tol 

382 # HACK: ECOS_BB uses ECOS' "tolerance for feasibility condition for 

383 # inaccurate solution" as the integrality tolerance. 

384 if self.ext.options.integrality_tol is not None: 

385 ecosOptions["feastol_inacc"] = self.ext.options.integrality_tol 

386 

387 # max_iterations 

388 if self.ext.options.max_iterations is not None: 

389 ecosOptions["max_iters"] = self.ext.options.max_iterations 

390 ecosOptions["mi_max_iters"] = self.ext.options.max_iterations 

391 

392 # Handle unsupported options. 

393 self._handle_unsupported_options( 

394 "lp_root_method", "lp_node_method", "timelimit", "treememory", 

395 "max_fsb_nodes", "hotstart") 

396 

397 # Assemble the solver input. 

398 arguments = {} 

399 arguments.update(self.int) 

400 arguments.update(ecosOptions) 

401 

402 # Debug print the solver input. 

403 if self._debug(): 

404 from pprint import pformat 

405 self._debug("Passing to ECOS:\n" + pformat(arguments)) 

406 

407 # Attempt to solve the problem. 

408 with self._header(), self._stopwatch(): 

409 solution = self.ecos.solve(**arguments) 

410 

411 # Debug print the solver output. 

412 if self._debug(): 

413 from pprint import pformat 

414 self._debug("Recevied from ECOS:\n" + pformat(solution)) 

415 

416 # Retrieve primals. 

417 primals = {} 

418 if self.ext.options.primals is not False: 

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

420 offset = self._ecosVarOffset[variable] 

421 dim = variable.dim 

422 value = list(solution["x"][offset:offset + dim]) 

423 primals[variable] = value 

424 

425 # Retrieve duals. 

426 duals = {} 

427 if self.ext.options.duals is not False: 

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

429 if isinstance(constraint, DummyConstraint): 

430 duals[constraint] = cvxopt.spmatrix( 

431 [], [], [], constraint.size) 

432 continue 

433 

434 dual = None 

435 indices = self._ecosConIndices[constraint] 

436 

437 if isinstance(constraint, AffineConstraint): 

438 if constraint.is_equality(): 

439 dual = list(-solution["y"][indices]) 

440 else: 

441 dual = list(solution["z"][indices]) 

442 elif isinstance(constraint, SOCConstraint): 

443 dual = solution["z"][indices] 

444 dual = list(dual) 

445 elif isinstance(constraint, RSOCConstraint): 

446 dual = solution["z"][indices] 

447 dual[1:] = -dual[1:] 

448 alpha = dual[0] - dual[-1] 

449 beta = dual[0] + dual[-1] 

450 z = list(-2.0 * dual[1:-1]) 

451 dual = [alpha, beta] + z 

452 elif isinstance(constraint, ExpConeConstraint): 

453 zxy = solution["z"][indices] 

454 dual = [zxy[1], zxy[2], zxy[0]] 

455 else: 

456 assert False, "Unexpected constraint type." 

457 

458 if type(dual) is list: 

459 if len(dual) == 1: 

460 dual = float(dual[0]) 

461 else: 

462 dual = cvxopt.matrix(dual, constraint.size) 

463 

464 duals[constraint] = dual 

465 

466 # Retrieve objective value. 

467 p = solution["info"]["pcost"] 

468 d = solution["info"]["dcost"] 

469 

470 if p is not None and d is not None: 

471 value = 0.5 * (p + d) 

472 elif p is not None: 

473 value = p 

474 elif d is not None: 

475 value = d 

476 else: 

477 value = None 

478 

479 if value is not None: 

480 # Add back the constant part. 

481 value += self.ext.no.function._constant_coef[0] 

482 

483 # Flip back the sign for maximization. 

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

485 value = -value 

486 

487 # Retrieve solution status. 

488 statusCode = solution["info"]["exitFlag"] 

489 if statusCode == 0: # optimal 

490 primalStatus = SS_OPTIMAL 

491 dualStatus = SS_OPTIMAL 

492 problemStatus = PS_FEASIBLE 

493 elif statusCode == 10: # inaccurate/optimal 

494 primalStatus = SS_OPTIMAL 

495 dualStatus = SS_OPTIMAL 

496 problemStatus = PS_FEASIBLE 

497 elif statusCode == 1: # primal infeasible 

498 primalStatus = SS_INFEASIBLE 

499 dualStatus = SS_UNKNOWN 

500 problemStatus = PS_INFEASIBLE 

501 elif statusCode == 11: # inaccurate/primal infeasible 

502 primalStatus = SS_INFEASIBLE 

503 dualStatus = SS_UNKNOWN 

504 problemStatus = PS_INFEASIBLE 

505 elif statusCode == 2: # dual infeasible 

506 primalStatus = SS_UNKNOWN 

507 dualStatus = SS_INFEASIBLE 

508 problemStatus = PS_UNBOUNDED 

509 elif statusCode == 12: # inaccurate/dual infeasible 

510 primalStatus = SS_UNKNOWN 

511 dualStatus = SS_INFEASIBLE 

512 problemStatus = PS_UNBOUNDED 

513 elif statusCode == -1: # iteration limit reached 

514 primalStatus = SS_PREMATURE 

515 dualStatus = SS_PREMATURE 

516 problemStatus = PS_UNKNOWN 

517 elif statusCode == -2: # search direction unreliable 

518 primalStatus = SS_UNKNOWN 

519 dualStatus = SS_UNKNOWN 

520 problemStatus = PS_UNSTABLE 

521 elif statusCode == -3: # numerically problematic 

522 primalStatus = SS_UNKNOWN 

523 dualStatus = SS_UNKNOWN 

524 problemStatus = PS_UNSTABLE 

525 elif statusCode == -4: # interrupted 

526 primalStatus = SS_PREMATURE 

527 dualStatus = SS_PREMATURE 

528 problemStatus = PS_UNKNOWN 

529 elif statusCode == -7: # solver error 

530 primalStatus = SS_FAILURE 

531 dualStatus = SS_FAILURE 

532 problemStatus = PS_UNKNOWN 

533 else: 

534 primalStatus = SS_UNKNOWN 

535 dualStatus = SS_UNKNOWN 

536 problemStatus = PS_UNKNOWN 

537 

538 return self._make_solution(value, primals, duals, primalStatus, 

539 dualStatus, problemStatus, {"ecos_sol": solution}) 

540 

541 

542# -------------------------------------- 

543__all__ = api_end(_API_START, globals())