Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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

2# Copyright (C) 2018-2019 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 

22 

23from ..apidoc import api_end, api_start 

24from ..constraints import (AffineConstraint, DummyConstraint, 

25 ExpConeConstraint, RSOCConstraint, SOCConstraint) 

26from ..expressions import AffineExpression, BinaryVariable, IntegerVariable 

27from ..modeling.footprint import Specification 

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

29 PS_UNKNOWN, PS_UNSTABLE, SS_FAILURE, 

30 SS_INFEASIBLE, SS_OPTIMAL, SS_PREMATURE, 

31 SS_UNKNOWN) 

32from .solver import Solver 

33 

34_API_START = api_start(globals()) 

35# ------------------------------- 

36 

37 

38class ECOSSolver(Solver): 

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

40 

41 SUPPORTED = Specification( 

42 objectives=[ 

43 AffineExpression], 

44 constraints=[ 

45 DummyConstraint, 

46 AffineConstraint, 

47 SOCConstraint, 

48 RSOCConstraint, 

49 ExpConeConstraint]) 

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 not in cls.SUPPORTED: 

59 if explain: 

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

61 else: 

62 return False 

63 

64 return (True, None) if explain else True 

65 

66 @classmethod 

67 def default_penalty(cls): 

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

69 return 1.0 # Stable free/open source solver. 

70 

71 @classmethod 

72 def test_availability(cls): 

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

74 cls.check_import("ecos") 

75 

76 @classmethod 

77 def names(cls): 

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

79 return "ecos", "ECOS", "Embedded Conic Solver" 

80 

81 @classmethod 

82 def is_free(cls): 

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

84 return True 

85 

86 def __init__(self, problem): 

87 """Initialize an ECOS solver interface. 

88 

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

90 """ 

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

92 

93 self._numVars = 0 

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

95 

96 self._ecosVarOffset = {} 

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

98 

99 self._ecosConIndices = dict() 

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

101 

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

103 

104 self._ecosModuleCache = None 

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

106 

107 def reset_problem(self): 

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

109 self.int = None 

110 

111 self._numVars = 0 

112 self._ecosVarOffset.clear() 

113 self._ecosConIndices.clear() 

114 

115 @property 

116 def ecos(self): 

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

118 

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

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

121 """ 

122 if self._ecosModuleCache is None: 

123 import ecos 

124 if hasattr(ecos, "ecos"): 

125 self._ecosModuleCache = ecos.ecos 

126 else: 

127 self._ecosModuleCache = ecos 

128 return self._ecosModuleCache 

129 

130 @property 

131 def array(self): 

132 """ECOS' array type.""" 

133 return self.ecos.np.array 

134 

135 @property 

136 def matrix(self): 

137 """ECOS' matrix type.""" 

138 return self.ecos.sparse.csc_matrix 

139 

140 def zeros(self, shape): 

141 """Create a zero array or a zero matrix, depending on ``shape``.""" 

142 if isinstance(shape, int) or len(shape) == 1: 

143 return self.ecos.np.zeros(shape) 

144 else: 

145 return self.matrix(shape) 

146 

147 def stack(self, *args): 

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

149 # In the case of matrices, stack vertically. 

150 if isinstance(args[0], self.ecos.sparse.base.spmatrix): 

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

152 assert isinstance(args[i], self.ecos.sparse.base.spmatrix) 

153 return self.ecos.sparse.vstack(args, format="csc") 

154 

155 # In the case of arrays, append them. 

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

157 assert isinstance(args[i], self.ecos.np.ndarray) 

158 return self.ecos.np.hstack(args) 

159 

160 def _affine_expression_to_G_and_h(self, expression): 

161 assert isinstance(expression, AffineExpression) 

162 

163 length = len(expression) 

164 

165 # Construct G. 

166 I, J, V = [], [], [] 

167 for variable, coefficients in expression._linear_coefs.items(): 

168 ecosVarOffset = self._ecosVarOffset[variable] 

169 

170 if not isinstance(coefficients, cvxopt.spmatrix): 

171 coefficients = cvxopt.sparse(coefficients) 

172 

173 I.extend(coefficients.I) 

174 J.extend([ecosVarOffset + j for j in coefficients.J]) 

175 V.extend(coefficients.V) 

176 G = self.matrix((V, (I, J)), (length, self._numVars)) 

177 

178 # Construct h. 

179 constant = expression._constant_coef 

180 if not isinstance(constant, cvxopt.matrix): 

181 constant = cvxopt.matrix(constant) 

182 h = self.array(constant, dtype=float).flatten() 

183 

184 return G, h 

185 

186 _Gh = _affine_expression_to_G_and_h 

187 

188 def _import_variables(self): 

189 offset = 0 

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

191 dim = variable.dim 

192 

193 # Mark integer variables. 

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

195 ecosIndices = range(offset, offset + dim) 

196 

197 if isinstance(variable, BinaryVariable): 

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

199 else: 

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

201 

202 # Register the variable. 

203 self._ecosVarOffset[variable] = offset 

204 offset += dim 

205 

206 # Import bounds. 

207 bounds = variable.bound_constraint 

208 if bounds: 

209 self._import_affine_constraint(bounds) 

210 

211 assert offset == self._numVars 

212 

213 def _append_equality(self, A, b): 

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

215 length = len(b) 

216 

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

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

219 

220 return range(offset, offset + length) 

221 

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

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

224 

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

226 if typecode == "l": 

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

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

229 elif typecode == "q": 

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

231 

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

233 length = len(h) 

234 

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

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

237 

238 if typecode == "q": 

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

240 elif typecode == "e": 

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

242 else: 

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

244 

245 return range(offset, offset + length) 

246 

247 def _import_affine_constraint(self, constraint): 

248 assert isinstance(constraint, AffineConstraint) 

249 

250 (G_smaller, h_smaller) = self._Gh(constraint.smaller) 

251 (G_greater, h_greater) = self._Gh(constraint.greater) 

252 G = G_smaller - G_greater 

253 h = h_greater - h_smaller 

254 

255 if constraint.is_equality(): 

256 return self._append_equality(G, h) 

257 else: 

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

259 

260 def _import_socone_constraint(self, constraint): 

261 assert isinstance(constraint, SOCConstraint) 

262 

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

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

265 

266 return self._append_inequality( 

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

268 

269 def _import_rscone_constraint(self, constraint): 

270 assert isinstance(constraint, RSOCConstraint) 

271 

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

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

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

275 

276 return self._append_inequality( 

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

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

279 

280 def _import_expcone_constraint(self, constraint): 

281 assert isinstance(constraint, ExpConeConstraint) 

282 

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

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

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

286 

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

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

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

290 # transform from our coordinates to theirs with the mapping 

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

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

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

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

295 return self._append_inequality( 

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

297 

298 def _import_objective(self): 

299 direction, objective = self.ext.no 

300 

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

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

303 

304 # Import coefficients. 

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

306 for localIndex in range(variable.dim): 

307 index = self._ecosVarOffset[variable] + localIndex 

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

309 

310 def _import_problem(self): 

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

312 

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

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

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

316 self.int = { 

317 # Objective function coefficients. 

318 "c": self.zeros(self._numVars), 

319 

320 # Linear equality left hand side. 

321 "A": self.matrix((0, self._numVars)), 

322 

323 # Linear equality right hand side. 

324 "b": self.array([]), 

325 

326 # Conic inequality left hand side. 

327 "G": self.matrix((0, self._numVars)), 

328 

329 # Conic inequality right hand side. 

330 "h": self.array([]), 

331 

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

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

334 

335 # Boolean variable indices. 

336 "bool_vars_idx": [], 

337 

338 # Integer variable indices. 

339 "int_vars_idx": [] 

340 } 

341 

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

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

344 

345 # Import variables with their bounds as affine constraints. 

346 self._import_variables() 

347 

348 # Import affine constraints. 

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

350 if isinstance(constraint, AffineConstraint): 

351 self._ecosConIndices[constraint] = \ 

352 self._import_affine_constraint(constraint) 

353 

354 # Import second order cone constraints. 

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

356 if isinstance(constraint, SOCConstraint): 

357 self._ecosConIndices[constraint] = \ 

358 self._import_socone_constraint(constraint) 

359 

360 # Import rotated second order cone constraints. 

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

362 if isinstance(constraint, RSOCConstraint): 

363 self._ecosConIndices[constraint] = \ 

364 self._import_rscone_constraint(constraint) 

365 

366 # Import exponential cone constraints. 

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

368 if isinstance(constraint, ExpConeConstraint): 

369 self._ecosConIndices[constraint] = \ 

370 self._import_expcone_constraint(constraint) 

371 

372 # Make sure that no unsupported constraints are present. 

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

374 assert isinstance(constraint, (DummyConstraint, AffineConstraint, 

375 SOCConstraint, RSOCConstraint, ExpConeConstraint)) 

376 

377 # Set objective. 

378 self._import_objective() 

379 

380 def _update_problem(self): 

381 raise NotImplementedError 

382 

383 def _solve(self): 

384 ecosOptions = {} 

385 

386 # verbosity 

387 beVerbose = self.verbosity() >= 1 

388 ecosOptions["verbose"] = beVerbose 

389 ecosOptions["mi_verbose"] = beVerbose 

390 

391 # rel_prim_fsb_tol, rel_dual_fsb_tol 

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

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

394 if feasibilityTols: 

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

396 

397 # abs_ipm_opt_tol 

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

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

400 

401 # rel_ipm_opt_tol 

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

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

404 

405 # abs_bnb_opt_tol 

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

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

408 

409 # rel_bnb_opt_tol 

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

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

412 

413 # integrality_tol 

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

415 # inaccurate solution" as the integrality tolerance. 

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

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

418 

419 # max_iterations 

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

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

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

423 

424 # Handle unsupported options. 

425 self._handle_unsupported_options( 

426 "lp_root_method", "lp_node_method", "timelimit", "treememory", 

427 "max_fsb_nodes", "hotstart") 

428 

429 # Assemble the solver input. 

430 arguments = {} 

431 arguments.update(self.int) 

432 arguments.update(ecosOptions) 

433 

434 # Debug print the solver input. 

435 if self._debug(): 

436 from pprint import pformat 

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

438 

439 # Attempt to solve the problem. 

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

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

442 

443 # Debug print the solver output. 

444 if self._debug(): 

445 from pprint import pformat 

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

447 

448 # Retrieve primals. 

449 primals = {} 

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

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

452 offset = self._ecosVarOffset[variable] 

453 dim = variable.dim 

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

455 primals[variable] = value 

456 

457 # Retrieve duals. 

458 duals = {} 

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

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

461 if isinstance(constraint, DummyConstraint): 

462 duals[constraint] = cvxopt.spmatrix( 

463 [], [], [], constraint.size) 

464 continue 

465 

466 dual = None 

467 indices = self._ecosConIndices[constraint] 

468 

469 if isinstance(constraint, AffineConstraint): 

470 if constraint.is_equality(): 

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

472 else: 

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

474 elif isinstance(constraint, SOCConstraint): 

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

476 dual = list(dual) 

477 elif isinstance(constraint, RSOCConstraint): 

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

479 dual[1:] = -dual[1:] 

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

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

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

483 dual = [alpha, beta] + z 

484 elif isinstance(constraint, ExpConeConstraint): 

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

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

487 else: 

488 assert False, "Unexpected constraint type." 

489 

490 if type(dual) is list: 

491 if len(dual) == 1: 

492 dual = float(dual[0]) 

493 else: 

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

495 

496 duals[constraint] = dual 

497 

498 # Retrieve objective value. 

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

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

501 

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

503 value = 0.5 * (p + d) 

504 elif p is not None: 

505 value = p 

506 elif d is not None: 

507 value = d 

508 else: 

509 value = None 

510 

511 if value is not None: 

512 # Add back the constant part. 

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

514 

515 # Flip back the sign for maximization. 

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

517 value = -value 

518 

519 # Retrieve solution status. 

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

521 if statusCode == 0: # optimal 

522 primalStatus = SS_OPTIMAL 

523 dualStatus = SS_OPTIMAL 

524 problemStatus = PS_FEASIBLE 

525 elif statusCode == 10: # inaccurate/optimal 

526 primalStatus = SS_OPTIMAL 

527 dualStatus = SS_OPTIMAL 

528 problemStatus = PS_FEASIBLE 

529 elif statusCode == 1: # primal infeasible 

530 primalStatus = SS_INFEASIBLE 

531 dualStatus = SS_UNKNOWN 

532 problemStatus = PS_INFEASIBLE 

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

534 primalStatus = SS_INFEASIBLE 

535 dualStatus = SS_UNKNOWN 

536 problemStatus = PS_INFEASIBLE 

537 elif statusCode == 2: # dual infeasible 

538 primalStatus = SS_UNKNOWN 

539 dualStatus = SS_INFEASIBLE 

540 problemStatus = PS_UNBOUNDED 

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

542 primalStatus = SS_UNKNOWN 

543 dualStatus = SS_INFEASIBLE 

544 problemStatus = PS_UNBOUNDED 

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

546 primalStatus = SS_PREMATURE 

547 dualStatus = SS_PREMATURE 

548 problemStatus = PS_UNKNOWN 

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

550 primalStatus = SS_UNKNOWN 

551 dualStatus = SS_UNKNOWN 

552 problemStatus = PS_UNSTABLE 

553 elif statusCode == -3: # numerically problematic 

554 primalStatus = SS_UNKNOWN 

555 dualStatus = SS_UNKNOWN 

556 problemStatus = PS_UNSTABLE 

557 elif statusCode == -4: # interrupted 

558 primalStatus = SS_PREMATURE 

559 dualStatus = SS_PREMATURE 

560 problemStatus = PS_UNKNOWN 

561 elif statusCode == -7: # solver error 

562 primalStatus = SS_FAILURE 

563 dualStatus = SS_FAILURE 

564 problemStatus = PS_UNKNOWN 

565 else: 

566 primalStatus = SS_UNKNOWN 

567 dualStatus = SS_UNKNOWN 

568 problemStatus = PS_UNKNOWN 

569 

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

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

572 

573 

574# -------------------------------------- 

575__all__ = api_end(_API_START, globals())