Coverage for picos/reforms/reform_dualize.py : 97.96%
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) 2020 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# ------------------------------------------------------------------------------
19"""Implementation of :class:`Dualization`."""
21from itertools import chain
23import cvxopt
25from ..apidoc import api_end, api_start
26from ..constraints import (AffineConstraint, DummyConstraint, LMIConstraint,
27 RSOCConstraint, SOCConstraint)
28from ..expressions import AffineExpression, Constant
29from ..expressions.algebra import rsoc, soc
30from ..expressions.variables import (CONTINUOUS_VARTYPES, RealVariable,
31 SymmetricVariable)
32from ..modeling import Problem
33from ..modeling.footprint import Footprint, Specification
34from ..modeling.solution import PS_INFEASIBLE, PS_UNBOUNDED, Solution
35from .reformulation import Reformulation
37_API_START = api_start(globals())
38# -------------------------------
41class Dualization(Reformulation):
42 """Lagrange dual problem reformulation."""
44 SUPPORTED = Specification(
45 objectives=[
46 AffineExpression],
47 variables=CONTINUOUS_VARTYPES,
48 constraints=[
49 DummyConstraint,
50 AffineConstraint,
51 SOCConstraint,
52 RSOCConstraint,
53 LMIConstraint])
55 @classmethod
56 def supports(cls, footprint):
57 """Implement :meth:`~.reformulation.Reformulation.supports`."""
58 return footprint.options.dualize and footprint in cls.SUPPORTED
60 @classmethod
61 def predict(cls, footprint):
62 """Implement :meth:`~.reformulation.Reformulation.predict`."""
63 vars = []
64 cons = []
66 # Objective direction.
67 dir = "min" if footprint.direction == "max" else "max"
69 # Objective function.
70 # HACK: Conversion adds a zero-fixed variable to the objective so that
71 # it is never constant. This makes prediction succeed in any case.
72 obj = AffineExpression.make_type(
73 shape=(1, 1), constant=False, nonneg=False)
75 # Variables and their bounds.
76 bound_cons = []
77 for vartype, k in footprint.variables.items():
78 cons.append((AffineConstraint.make_type(
79 dim=vartype.subtype.dim, eq=True), k))
81 if vartype.subtype.bnd:
82 bound_cons.append((AffineConstraint.make_type(
83 dim=vartype.subtype.bnd, eq=False), k))
85 # Constraints by type.
86 for contype, k in chain(footprint.constraints.items(), bound_cons):
87 if contype.clstype is DummyConstraint:
88 pass
89 elif contype.clstype is AffineConstraint:
90 d = contype.subtype.dim
92 # Dual variable.
93 vars.append((RealVariable.make_var_type(dim=d, bnd=0), k))
95 # Conic membership constraint.
96 if not contype.subtype.eq:
97 cons.append((
98 AffineConstraint.make_type(dim=d, eq=False), k))
99 elif contype.clstype is SOCConstraint:
100 n = contype.subtype.argdim
101 d = n + 1
103 # Dual variable.
104 vars.append((RealVariable.make_var_type(dim=d, bnd=0), k))
106 # Conic membership constraint.
107 cons.append((SOCConstraint.make_type(argdim=n), k))
108 elif contype.clstype is RSOCConstraint:
109 n = contype.subtype.argdim
110 d = n + 2
112 # Dual variable.
113 vars.append((RealVariable.make_var_type(dim=d, bnd=0), k))
115 # Conic membership constraint.
116 cons.append((RSOCConstraint.make_type(argdim=n), k))
117 elif contype.clstype is LMIConstraint:
118 n = contype.subtype.diag
119 d = n*(n + 1)//2
121 # Dual variable.
122 vars.append((SymmetricVariable.make_var_type(dim=d, bnd=0), k))
124 # Conic membership constraint.
125 cons.append((LMIConstraint.make_type(diag=n), k))
126 else:
127 assert False, "Unexpected constraint type."
129 # HACK: See above.
130 vars.append((RealVariable.make_var_type(dim=1, bnd=2), 1))
132 # Option changes.
133 nd_opts = footprint.options.updated(
134 dualize=False,
135 primals=footprint.options.duals,
136 duals=footprint.options.primals
137 ).nondefaults
139 return Footprint.from_types(dir, obj, vars, cons, nd_opts)
141 def forward(self):
142 """Implement :meth:`~.reformulation.Reformulation.forward`."""
143 self._pc2dv = {} # Maps primal constraints to dual variables.
144 self._pv2dc = {} # Maps primal variables to dual constraints.
146 P = self.input
148 # Create the dual problem from scratch.
149 D = self.output = Problem(copyOptions=P.options)
151 # Reset the 'dualize' option so that solvers can load the dual.
152 D.options.dualize = False
154 # Swap 'primals' and 'duals' option.
155 D.options.primals = P.options.duals
156 D.options.duals = P.options.primals
158 # Retrieve the primal objective in standard form.
159 original_primal_direction, primal_objective = P.objective.normalized
160 if original_primal_direction == "max":
161 primal_objective = -primal_objective
163 # Start building the dual objective function.
164 dual_objective = primal_objective.cst
166 # Start building the dual linear equalities for each primal variable.
167 obj_coefs = primal_objective._linear_coefs
168 dual_equality_terms = {var: -Constant(obj_coefs[var].T)
169 if var in obj_coefs else AffineExpression.zero(var.dim)
170 for var in P.variables.values()}
172 # Turn variable bounds into affine inequality constraints.
173 # NOTE: If bound-free variables are required by another reformulation or
174 # solver in the future, there should be a reformulation for this.
175 bound_cons = []
176 for primal_var in P.variables.values():
177 bound_constraint = primal_var.bound_constraint
178 if bound_constraint:
179 bound_cons.append(bound_constraint)
181 # Handle the primal constraints.
182 # NOTE: Constraint conversion works as follows:
183 # 1. Transform constraint into conic standard form:
184 # fx = Ax-b >_K 0.
185 # 2. Create a dual variable for the constraint.
186 # 3. Constrain the dual variable to be a member of the dual
187 # cone K*.
188 # 4. Extend the dual's linear equalities for each primal variable.
189 # 5. Extend the dual's objective function.
190 # TODO: Consider adding a ConicConstraint abstract subclass of
191 # Constraint with a method conic_form that returns a pair of cone
192 # member (affine expression) and cone. The conic standard form
193 # function f(x) and the dual cone can then be obtained by negation
194 # and through a new Cone.dual method, respectively.
195 for primal_con in chain(P.constraints.values(), bound_cons):
196 if isinstance(primal_con, DummyConstraint):
197 pass
198 elif isinstance(primal_con, AffineConstraint):
199 # Transform constraint into conic standard form.
200 fx = primal_con.ge0.vec
202 # Create a dual variable for the constraint.
203 self._pc2dv[primal_con] = dual_var = RealVariable(
204 "aff_{}".format(primal_con.id), len(fx))
206 # Constrain the dual variable to be a member of the dual cone.
207 # NOTE: In the equality case, the dual cone is the real field.
208 if primal_con.is_inequality():
209 D.add_constraint(dual_var >= 0)
211 # Extend the dual's linear equalities for each primal variable.
212 for var, coef in fx._linear_coefs.items():
213 dual_equality_terms[var] += coef.T*dual_var
215 # Extend the dual's objective function.
216 dual_objective -= fx._constant_coef.T*dual_var
217 elif isinstance(primal_con, SOCConstraint):
218 # Transform constraint into conic standard form.
219 fx = (primal_con.ub // primal_con.ne.vec)
221 # Create a dual variable for the constraint.
222 self._pc2dv[primal_con] = dual_var = RealVariable(
223 "soc_{}".format(primal_con.id), len(fx))
225 # Constrain the dual variable to be a member of the dual cone.
226 D.add_constraint(dual_var << soc())
228 # Extend the dual's linear equalities for each primal variable.
229 for var, coef in fx._linear_coefs.items():
230 dual_equality_terms[var] += coef.T*dual_var
232 # Extend the dual's objective function.
233 dual_objective -= fx._constant_coef.T*dual_var
234 elif isinstance(primal_con, RSOCConstraint):
235 # Transform constraint into conic standard form.
236 fx = (primal_con.ub1 // primal_con.ub2 // primal_con.ne.vec)
238 # Create a dual variable for the constraint.
239 self._pc2dv[primal_con] = dual_var = RealVariable(
240 "rso_{}".format(primal_con.id), len(fx))
242 # Constrain the dual variable to be a member of the dual cone.
243 D.add_constraint(dual_var << rsoc(4))
245 # Extend the dual's linear equalities for each primal variable.
246 for var, coef in fx._linear_coefs.items():
247 dual_equality_terms[var] += coef.T*dual_var
249 # Extend the dual's objective function.
250 dual_objective -= fx._constant_coef.T*dual_var
251 elif isinstance(primal_con, LMIConstraint):
252 # Transform constraint into conic standard form.
253 fx = primal_con.psd
255 # Create a dual variable for the constraint.
256 self._pc2dv[primal_con] = dual_var = SymmetricVariable(
257 "lmi_{}".format(primal_con.id), fx.shape)
259 # Constrain the dual variable to be a member of the dual cone.
260 D.add_constraint(dual_var >> 0)
262 # Extend the dual's linear equalities for each primal variable.
263 for var, coef in fx._linear_coefs.items():
264 dual_equality_terms[var] += coef.T*dual_var.vec
266 # Extend the dual's objective function.
267 dual_objective -= fx._constant_coef.T*dual_var.vec
268 else:
269 assert False, "Unexpected constraint type."
271 # Add the finished linear equalities for each variable.
272 for var, term in dual_equality_terms.items():
273 self._pv2dc[var] = D.add_constraint(term == 0)
275 # Adjust dual optimization direction to obtain the duality property.
276 if original_primal_direction == "max":
277 dual_direction = "min"
278 dual_objective = -dual_objective
279 else:
280 dual_direction = "max"
282 # HACK: Add a zero-fixed variable to the objective so that it is never
283 # constant. This makes prediction succeed in any case.
284 dual_objective += RealVariable("zero", 1, lower=0, upper=0)
286 # Set the finished dual objective.
287 D.objective = dual_direction, dual_objective
289 def update(self):
290 """Implement :meth:`~.reformulation.Reformulation.update`."""
291 # TODO: Implement updates to an existing dualization. This is optional.
292 raise NotImplementedError
294 def backward(self, solution):
295 """Implement :meth:`~.reformulation.Reformulation.backward`."""
296 P = self.input
298 # Require primal values to be properly vectorized.
299 if not solution.vectorizedPrimals:
300 raise NotImplementedError(
301 "Solving with dualize=True is not supported with solvers that "
302 "produce unvectorized primal solutions.")
304 # Retrieve primals for the primal problem.
305 primals = {}
306 for var in P.variables.values():
307 con = self._pv2dc[var]
309 if con in solution.duals and solution.duals[con] is not None:
310 primal = -solution.duals[con]
311 else:
312 primal = None
314 if primal is not None:
315 primals[var] = cvxopt.matrix(primal)
316 else:
317 primals[var] = None
319 # Retrieve duals for the primal problem.
320 duals = {}
321 for con in P.constraints.values():
322 if isinstance(con, DummyConstraint):
323 continue
325 var = self._pc2dv[con]
326 dual = solution.primals[var] if var in solution.primals else None
328 if dual is not None:
329 duals[con] = var._vec.devectorize(cvxopt.matrix(dual))
330 else:
331 duals[con] = None
333 # Swap infeasible/unbounded for a problem status.
334 if solution.problemStatus == PS_UNBOUNDED:
335 problemStatus = PS_INFEASIBLE
336 elif solution.problemStatus == PS_INFEASIBLE:
337 problemStatus = PS_UNBOUNDED
338 else:
339 problemStatus = solution.problemStatus
341 return Solution(
342 primals=primals,
343 duals=duals,
344 problem=solution.problem,
345 solver=solution.solver,
346 primalStatus=solution.dualStatus, # Intentional swap.
347 dualStatus=solution.primalStatus, # Intentional swap.
348 problemStatus=problemStatus,
349 searchTime=solution.searchTime,
350 info=solution.info,
351 vectorizedPrimals=True,
352 reportedValue=solution.reportedValue)
355# --------------------------------------
356__all__ = api_end(_API_START, globals())