Coverage for picos/solvers/solver

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

5# This file is part of PICOS.

7# PICOS is free software: you can redistribute it and/or modify it under the

8# terms of the GNU General Public License as published by the Free Software

9# Foundation, either version 3 of the License, or (at your option) any later

10# version.

11#

12# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY

13# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR

14# A PARTICULAR PURPOSE. See the GNU General Public License for more details.

15#

16# You should have received a copy of the GNU General Public License along with

17# this program. If not, see <http://www.gnu.org/licenses/>.

18# ------------------------------------------------------------------------------

20"""Implementation of :class:`CVXOPTSolver`."""

22import cvxopt

23import numpy

25from ..apidoc import api_end, api_start

26from ..constraints import (AffineConstraint, DummyConstraint, LMIConstraint,

27 LogSumExpConstraint, RSOCConstraint, SOCConstraint)

28from ..expressions import CONTINUOUS_VARTYPES, AffineExpression, LogSumExp

29from ..modeling.footprint import Specification

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

31 PS_UNKNOWN, SS_INFEASIBLE, SS_OPTIMAL,

32 SS_UNKNOWN, Solution)

33from .solver import Solver

35_API_START = api_start(globals())

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

39class CVXOPTSolver(Solver):

40 """Interface to the CVXOPT solver.

42 Also used as an interface to the

43 :class:`SMCP solver <picos.solvers.solver_smcp.SMCPSolver>`.

44 """

46 SUPPORTED = Specification(

47 objectives=[

48 AffineExpression,

49 LogSumExp],

50 variables=CONTINUOUS_VARTYPES,

51 constraints=[

52 DummyConstraint,

53 AffineConstraint,

54 SOCConstraint,

55 RSOCConstraint,

56 LMIConstraint,

57 LogSumExpConstraint])

59 @classmethod

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

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

62 result = Solver.supports(footprint, explain)

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

64 return result

66 if footprint not in cls.SUPPORTED:

67 if explain:

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

69 else:

70 return False

72 return (True, None) if explain else True

74 @classmethod

75 def default_penalty(cls):

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

77 return 1.0 # Stable free/open source solver.

79 @classmethod

80 def test_availability(cls):

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

82 # CVXOPT is a dependency of PICOS, so it is always available.

83 pass

85 @classmethod

86 def names(cls):

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

88 return "cvxopt", "CVXOPT", "Python Convex Optimization Solver", None

90 @classmethod

91 def is_free(cls):

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

93 return True

95 @property

96 def is_smcp(self):

97 """Whether to implement SMCP instead of CVXOPT."""

98 return False

100 def __init__(self, problem):

101 """Initialize a CVXOPT solver interface.

102

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

104 """

105 super(CVXOPTSolver, self).__init__(problem)

106

107 self._numVars = 0

108 """Total number of scalar variables passed to CVXOPT."""

109

110 self._cvxoptVarOffset = {}

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

112

113 # HACK: Setting this to false would result in variable bounds to be

114 # ignored, instead of added to the linear inequalities matrix.

115 # This is used by a Problem method that prints the problem to a

116 # file using a CVXOPTSolver instance.

117 self.import_variable_bounds = True

118

119 # SMCP currently offers no option reset, so keep a backup.

120 if self.is_smcp:

121 import smcp

122 self._smcp_default_options = smcp.solvers.options.copy()

123

124 def reset_problem(self):

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

126 self.int = None

127

128 self._numVars = 0

129 self._cvxoptVarOffset.clear()

130

131 def _is_geometric_program(self):

132 return isinstance(self.ext.no.function, LogSumExp) or \

133 any(isinstance(constraint, LogSumExpConstraint)

134 for constraint in self.ext.constraints.values())

135

136 def _affine_expression_to_G_and_h(self, expression):

137 assert isinstance(expression, AffineExpression)

138

139 return expression.sparse_matrix_form(

140 varOffsetMap=self._cvxoptVarOffset, dense_b=True)

141

142 _Gh = _affine_expression_to_G_and_h

143

144 def _import_variables_without_bounds(self):

145 offset = 0

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

147 dim = variable.dim

148

149 # Register the variable.

150 self._cvxoptVarOffset[variable] = offset

151 offset += dim

152

153 assert offset == self._numVars

154

155 def _import_variable_bounds(self):

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

157 if self.import_variable_bounds:

158 bounds = variable.bound_constraint

159 if bounds:

160 self._import_affine_constraint(bounds)

161

162 def _import_objective(self):

163 direction, objective = self.ext.no

164

165 if self._is_geometric_program():

166 if isinstance(objective, LogSumExp):

167 (F, g) = self._Gh(objective.x)

168 else:

169 assert isinstance(objective, AffineExpression)

170 (F, g) = self._Gh(objective)

171

172 # NOTE: Needs to happen before LogSumExpConstraint are added.

173 self.int["K"] = [F.size[0]]

174

175 if direction == "max":

176 F, g = -F, -g

177

178 self.int["F"] = F

179 self.int["g"] = g

180 else:

181 (c, _) = self._Gh(objective)

182

183 if direction == "max":

184 c = -c

185

186 # Must be a dense column-vector.

187 c = cvxopt.matrix(c).T

188

189 self.int["c"] = c

190

191 def _import_affine_constraint(self, constraint):

192 assert isinstance(constraint, AffineConstraint)

193

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

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

196

197 G = G_smaller - G_greater

198 h = h_greater - h_smaller

199

200 if constraint.is_equality():

201 self.int["A"] = cvxopt.sparse([self.int["A"], G])

202 self.int["b"] = cvxopt.matrix([self.int["b"], h])

203 else:

204 self.int["Gl"] = cvxopt.sparse([self.int["Gl"], G])

205 self.int["hl"] = cvxopt.matrix([self.int["hl"], h])

206

207 def _import_soc_constraint(self, constraint):

208 assert isinstance(constraint, SOCConstraint)

209

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

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

212

213 self.int["Gq"].append(cvxopt.sparse([-c, -A]))

214 self.int["hq"].append(cvxopt.matrix([d, b]))

215

216 def _import_rsoc_constraint(self, constraint):

217 assert isinstance(constraint, RSOCConstraint)

218

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

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

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

222

223 self.int["Gq"].append(cvxopt.sparse([-c1 - c2, -2 * A, c2 - c1]))

224 self.int["hq"].append(cvxopt.matrix([d1 + d2, 2 * b, d1 - d2]))

225

226 def _import_lmi_constraint(self, constraint):

227 assert isinstance(constraint, LMIConstraint)

228

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

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

231

232 self.int["Gs"].append(G_smaller - G_greater)

233 self.int["hs"].append(h_greater - h_smaller)

234

235 def _import_lse_constraint(self, constraint):

236 assert isinstance(constraint, LogSumExpConstraint)

237

238 (F, g) = self._Gh(constraint.le0.x)

239

240 self.int["F"] = cvxopt.sparse([self.int["F"], F])

241 self.int["g"] = cvxopt.matrix([self.int["g"], g])

242 self.int["K"].append(F.size[0])

243

244 def _import_constraint(self, constraint):

245 if isinstance(constraint, AffineConstraint):

246 self._import_affine_constraint(constraint)

247 elif isinstance(constraint, SOCConstraint):

248 self._import_soc_constraint(constraint)

249 elif isinstance(constraint, RSOCConstraint):

250 self._import_rsoc_constraint(constraint)

251 elif isinstance(constraint, LMIConstraint):

252 self._import_lmi_constraint(constraint)

253 elif isinstance(constraint, LogSumExpConstraint):

254 self._import_lse_constraint(constraint)

255 else:

256 assert isinstance(constraint, DummyConstraint), \

257 "Unexpected constraint type: {}".format(

258 constraint.__class__.__name__)

259

260 def _import_problem(self):

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

262

263 # CVXOPT's internal problem representation is stateless; a number of

264 # matrices are supplied to the appropriate solver function each time a

265 # search is started. These matrices are thus stored in self.int.

266 self.int = {

267 # Objective function coefficients.

268 "c": None,

269

270 # Linear equality left hand side.

271 "A": cvxopt.spmatrix([], [], [], (0, self._numVars), tc="d"),

272

273 # Linear equality right hand side.

274 "b": cvxopt.matrix([], (0, 1), tc="d"),

275

276 # Linear inequality left hand side.

277 "Gl": cvxopt.spmatrix([], [], [], (0, self._numVars), tc="d"),

278

279 # Linear inequality right hand side.

280 "hl": cvxopt.matrix([], (0, 1), tc="d"),

281

282 # Second order cone inequalities left hand sides.

283 "Gq": [],

284

285 # Second order cone inequalities right hand sides.

286 "hq": [],

287

288 # Semidefinite cone inequalities left hand sides.

289 "Gs": [],

290

291 # Semidefinite cone inequalities right hand sides.

292 "hs": [],

293

294 # Geometric program data.

295 "F": None,

296 "g": None,

297 "K": None

298 }

299

300 # Import variables without their bounds.

301 self._import_variables_without_bounds()

302

303 # Set objective.

304 # NOTE: This needs to happen before constraints are added as

305 # self.int["K"][0] refers to an LSE objective while

306 # self.int["K"][i] with i > 0 refers to LSE constraints.

307 self._import_objective()

308

309 # Import constraints.

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

311 self._import_constraint(constraint)

312

313 # Import variable bounds as additional affine constraints.

314 # NOTE: This needs to happen after constraints are added due to how the

315 # dual values for constraints are extracted.

316 self._import_variable_bounds()

317

318 def _update_problem(self):

319 raise NotImplementedError

320

321 def _solve(self):

322 if self.is_smcp:

323 import smcp

324

325 p = self.int

326 isGP = self._is_geometric_program()

327

328 # Clear all options set previously. This is necessary because CVXOPT

329 # options are global, and might be changed even by another problem.

330 cvxopt.solvers.options.clear()

331

332 # verbosity

333 cvxopt.solvers.options["show_progress"] = (self.verbosity() >= 1)

334

335 # rel_prim_fsb_tol, rel_dual_fsb_tol

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

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

338 if feasibilityTols:

339 cvxopt.solvers.options["feastol"] = min(feasibilityTols)

340

341 # abs_ipm_opt_tol

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

343 cvxopt.solvers.options["abstol"] = self.ext.options.abs_ipm_opt_tol

344

345 # rel_ipm_opt_tol

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

347 cvxopt.solvers.options["reltol"] = self.ext.options.rel_ipm_opt_tol

348

349 # max_iterations

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

351 cvxopt.solvers.options["maxiters"] = self.ext.options.max_iterations

352 else:

353 cvxopt.solvers.options["maxiters"] = int(1e6)

354

355 # cvxopt_kktsolver

356 if self.ext.options.cvxopt_kktsolver is not None:

357 userKKT = self.ext.options.cvxopt_kktsolver

358 else:

359 userKKT = None

360

361 # cvxopt_kktreg

362 if self.ext.options.cvxopt_kktreg is not None:

363 cvxopt.solvers.options["kktreg"] = self.ext.options.cvxopt_kktreg

364

365 # Handle unsupported options.

366 self._handle_unsupported_options(

367 "lp_root_method", "lp_node_method", "timelimit", "treememory",

368 "max_fsb_nodes", "hotstart")

369

370 # TODO: Add CVXOPT-sepcific options. Candidates are:

371 # - refinement

372

373 # Set options for SMCP.

374 if self.is_smcp:

375 # Restore default options.

376 smcp.solvers.options = self._smcp_default_options.copy()

377

378 # Copy options also used by CVXOPT.

379 smcp.solvers.options.update(cvxopt.solvers.options)

380

381 # Further handle "verbose" option.

382 smcp.solvers.options["debug"] = (self.verbosity() >= 2)

383

384 # TODO: Add SMCP-sepcific options.

385

386 if self._debug():

387 from pprint import pformat

388 self._debug("Setting options:\n{}\n".format(pformat(

389 smcp.solvers.options if self.is_smcp

390 else cvxopt.solvers.options)))

391

392 # Print a header.

393 if self.is_smcp:

394 subsolverText = None

395 else:

396 if isGP:

397 subsolverText = "internal GP solver"

398 else:

399 subsolverText = "internal CONELP solver"

400

401 # Further prepare the problem for the CVXOPT/SMCP CONELP solvers.

402 # TODO: This should be done during import.

403 if not isGP:

404 # Retrieve the structure of the cone, which is a cartesian product

405 # of the non-negative orthant of dimension l, a number of second

406 # order cones with dimensions in q and a number of positive

407 # semidefinite cones with dimensions in s.

408 dims = {

409 "l": p["Gl"].size[0],

410 "q": [Gqi.size[0] for Gqi in p["Gq"]],

411 "s": [int(numpy.sqrt(Gsi.size[0])) for Gsi in p["Gs"]]

412 }

413

414 # Construct G and h to contain all conic inequalities, starting with

415 # those with respect to the non-negative orthant.

416 G = p["Gl"]

417 h = p["hl"]

418

419 # SMCP's ConeLP solver does not handle (linear) equalities, so cast

420 # them as inequalities.

421 if self.is_smcp:

422 smcp_eps = min(feasibilityTols)

423 if p["A"].size[0] > 0:

424 G = cvxopt.sparse([G, p["A"]])

425 G = cvxopt.sparse([G, -p["A"]])

426 h = cvxopt.matrix([h, p["b"]+smcp_eps])

427 h = cvxopt.matrix([h, smcp_eps-p["b"]])

428 dims["l"] += (2 * p["A"].size[0])

429

430 # Remove the lines in G and h corresponding to 0==0 or 0<=0

431 JP = list(set(G.I))

432 IP = range(len(JP))

433 VP = [1] * len(JP)

434

435 if len(JP) != dims["l"]:

436 # is there a constraint of the form 0<=a, (a<0) ?

437 if any([b < -smcp_eps for (i, b) in enumerate(h)

438 if i not in JP]):

439 raise Exception(

440 'infeasible constraint of the form '

441 '0 <= a, with a<0')

442

443 # left-multiply with PPP-matrix to remove 0-constraints

444 PPP = cvxopt.spmatrix(VP, IP, JP, (len(IP), G.size[0]))

445 dims["l"] = len(JP)

446 G = PPP * G

447 h = PPP * h

448

449 # Add second-order cone inequalities.

450 for i in range(len(dims["q"])):

451 G = cvxopt.sparse([G, p["Gq"][i]])

452 h = cvxopt.matrix([h, p["hq"][i]])

453

454 # Add semidefinite cone inequalities.

455 for i in range(len(dims["s"])):

456 G = cvxopt.sparse([G, p["Gs"][i]])

457 h = cvxopt.matrix([h, p["hs"][i]])

458

459 # Remove zero lines from linear equality constraint matrix, as

460 # CVXOPT expects this matrix to have full row rank.

461 JP = list(set(p["A"].I))

462 IP = range(len(JP))

463 VP = [1] * len(JP)

464

465 # Skip solution on an infeasible constraint.

466 if any([b for (i, b) in enumerate(p["b"]) if i not in JP]):

467 return Solution(

468 primals=None, duals=None, problem=self.ext, solver="PICOS",

469 primalStatus=SS_INFEASIBLE, dualStatus=SS_UNKNOWN,

470 problemStatus=PS_INFEASIBLE, vectorizedPrimals=True)

471

472 P = cvxopt.spmatrix(VP, IP, JP, (len(IP), p["A"].size[0]))

473 A = P * p["A"]

474 b = P * p["b"]

475

476 # Attempt to solve the problem.

477 with self._header(subsolverText), self._stopwatch():

478 if self.is_smcp:

479 if self._debug():

480 self._debug("Calling smcp.solvers.conelp(c, G, h, dims) "

481 "with\nc:\n{}\nG:\n{}\nh:\n{}\ndims:\n{}\n"

482 .format(p["c"], G, h, dims))

483 try:

484 result = smcp.solvers.conelp(p["c"], G, h, dims)

485 except TypeError:

486 # HACK: Work around "'NoneType' object is not subscriptable"

487 # exception with infeasible/unbounded problems.

488 result = None

489 else:

490 kwargs = {}

491

492 if userKKT:

493 kwargs["kktsolver"] = userKKT

494 elif isGP:

495 # Use the more reliable LDL solver right away.

496 kwargs["kktsolver"] = "ldl"

497 else:

498 # Try the fast but unreliable CHOL solver first.

499 kwargs["kktsolver"] = "chol"

500

501 if isGP:

502 result = cvxopt.solvers.gp(p["K"], p["F"], p["g"], p["Gl"],

503 p["hl"], p["A"], p["b"], **kwargs)

504 else:

505 try:

506 result = cvxopt.solvers.conelp(

507 p["c"], G, h, dims, A, b, **kwargs)

508

509 if not userKKT and result["status"] == "unknown":

510 raise ValueError("The first solution attempt with "

511 "CHOL as a KKT solver returnd a solution with "

512 "unknown status. This exception triggers "

513 "another solution attempt using LDL.")

514 except ValueError:

515 # NOTE: Apart from the one created by PICOS above, there

516 # are at least two more ValueError produced by

517 # CVXOPT when an unreliable KKT solver is used.

518 if userKKT:

519 raise # Always respect the user's choice.

520 else:

521 # Re-solve using the LDL solver.

522 # TODO: Consider pre-solving on PICOS end to prevent

523 # the "Rank(A) < p or Rank([G; A]) < n" error.

524 kwargs["kktsolver"] = "ldl"

525 result = cvxopt.solvers.conelp(

526 p["c"], G, h, dims, A, b, **kwargs)

527

528 # Retrieve primals.

529 primals = {}

530 if self.ext.options.primals is not False and result is not None \

531 and result["x"] is not None:

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

533 offset = self._cvxoptVarOffset[variable]

534 value = list(result["x"][offset:offset + variable.dim])

535 primals[variable] = value

536

537 # Retrieve duals.

538 duals = {}

539 if self.ext.options.duals is not False and result is not None:

540 (indy, indzl, indzq, indznl, indzs) = (0, 0, 0, 0, 0)

541

542 if isGP:

543 zkey = "zl"

544 zqkey = "zq"

545 zskey = "zs"

546 else:

547 zkey = "z"

548 zqkey = "z"

549 zskey = "z"

550 indzq = dims["l"]

551 indzs = dims["l"] + sum(dims["q"])

552

553 if self.is_smcp:

554 # Equality constraints were cast as two inequalities.

555 ieq = p["Gl"].size[0]

556 neq = (dims["l"] - ieq) // 2

557 soleq = result["z"][ieq:ieq + neq]

558 soleq -= result["z"][ieq + neq:ieq + 2 * neq]

559 else:

560 soleq = result["y"]

561

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

563 if isinstance(constraint, DummyConstraint):

564 duals[constraint] = cvxopt.spmatrix(

565 [], [], [], constraint.size)

566 continue

567

568 dual = None

569 consSz = len(constraint)

570

571 if isinstance(constraint, AffineConstraint):

572 if constraint.is_equality():

573 if soleq is not None:

574 dual = -(P.T * soleq)[indy:indy + consSz]

575 indy += consSz

576 else:

577 if result[zkey] is not None:

578 dual = result[zkey][indzl:indzl + consSz]

579 indzl += consSz

580 elif isinstance(constraint, SOCConstraint) \

581 or isinstance(constraint, RSOCConstraint):

582 if result[zqkey] is not None:

583 if isGP:

584 dual = result[zqkey][indzq]

585 dual[1:] = -dual[1:]

586 indzq += 1

587 else:

588 dual = result[zqkey][indzq:indzq + consSz]

589 if isinstance(constraint, RSOCConstraint):

590 # RScone were cast as a SOcone on import, so

591 # transform the dual to a proper RScone dual.

592 alpha = dual[0] + dual[-1]

593 beta = dual[0] - dual[-1]

594 z = 2.0 * dual[1:-1]

595 dual = cvxopt.matrix([alpha, beta, z])

596 indzq += consSz

597 elif isinstance(constraint, LMIConstraint):

598 if result[zskey] is not None:

599 matSz = constraint.size[0]

600 if isGP:

601 dual = cvxopt.matrix(

602 result[zskey][indzs], (matSz, matSz))

603 indzs += 1

604 else:

605 dual = cvxopt.matrix(

606 result[zskey][indzs:indzs + consSz],

607 (matSz, matSz))

608 indzs += consSz

609 elif isinstance(constraint, LogSumExpConstraint):

610 # TODO: Retrieve LSE duals.

611 indznl += 1

612

613 duals[constraint] = dual

614

615 # Retrieve objective value.

616 if result is None:

617 value = None

618 elif isGP:

619 value = None

620 else:

621 p = result['primal objective']

622 d = result['dual objective']

623

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

625 value = 0.5 * (p + d)

626 elif p is not None:

627 value = p

628 elif d is not None:

629 value = d

630 else:

631 value = None

632

633 if value is not None and self.ext.no.direction == "max":

634 value = -value

635

636 if self.is_smcp:

637 value = -value

638

639 # Retrieve solution status.

640 status = result["status"] if result else "unknown"

641 if status == "optimal":

642 primalStatus = SS_OPTIMAL

643 dualStatus = SS_OPTIMAL

644 problemStatus = PS_FEASIBLE

645 elif status == "primal infeasible":

646 primalStatus = SS_INFEASIBLE

647 dualStatus = SS_UNKNOWN

648 problemStatus = PS_INFEASIBLE

649 elif status == "dual infeasible":

650 primalStatus = SS_UNKNOWN

651 dualStatus = SS_INFEASIBLE

652 problemStatus = PS_UNBOUNDED

653 else:

654 primalStatus = SS_UNKNOWN

655 dualStatus = SS_UNKNOWN

656 problemStatus = PS_UNKNOWN

657

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

659 dualStatus, problemStatus, {"cvxopt_sol": result})

660

661

662# --------------------------------------

663__all__ = api_end(_API_START, globals())

Coverage for picos/solvers/solver_cvxopt.py: 83.14%

338 statements