Coverage for picos/reforms/reform_dualize.py : 97.87%
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, LMIConstraint, RSOCConstraint,
27 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 AffineConstraint,
50 SOCConstraint,
51 RSOCConstraint,
52 LMIConstraint])
54 @classmethod
55 def supports(cls, footprint):
56 """Implement :meth:`~.reformulation.Reformulation.supports`."""
57 return footprint.options.dualize and footprint in cls.SUPPORTED
59 @classmethod
60 def predict(cls, footprint):
61 """Implement :meth:`~.reformulation.Reformulation.predict`."""
62 vars = []
63 cons = []
65 # Objective direction.
66 dir = "min" if footprint.direction == "max" else "max"
68 # Objective function.
69 # HACK: Conversion adds a zero-fixed variable to the objective so that
70 # it is never constant. This makes prediction succeed in any case.
71 obj = AffineExpression.make_type(
72 shape=(1, 1), constant=False, nonneg=False)
74 # Variables and their bounds.
75 bound_cons = []
76 for vartype, k in footprint.variables.items():
77 cons.append((AffineConstraint.make_type(
78 dim=vartype.subtype.dim, eq=True), k))
80 if vartype.subtype.bnd:
81 bound_cons.append((AffineConstraint.make_type(
82 dim=vartype.subtype.bnd, eq=False), k))
84 # Constraints by type.
85 for contype, k in chain(footprint.constraints.items(), bound_cons):
86 if contype.clstype is AffineConstraint:
87 d = contype.subtype.dim
89 # Dual variable.
90 vars.append((RealVariable.make_var_type(dim=d, bnd=0), k))
92 # Conic membership constraint.
93 if not contype.subtype.eq:
94 cons.append((
95 AffineConstraint.make_type(dim=d, eq=False), k))
96 elif contype.clstype is SOCConstraint:
97 n = contype.subtype.argdim
98 d = n + 1
100 # Dual variable.
101 vars.append((RealVariable.make_var_type(dim=d, bnd=0), k))
103 # Conic membership constraint.
104 cons.append((SOCConstraint.make_type(argdim=n), k))
105 elif contype.clstype is RSOCConstraint:
106 n = contype.subtype.argdim
107 d = n + 2
109 # Dual variable.
110 vars.append((RealVariable.make_var_type(dim=d, bnd=0), k))
112 # Conic membership constraint.
113 cons.append((RSOCConstraint.make_type(argdim=n), k))
114 elif contype.clstype is LMIConstraint:
115 n = contype.subtype.diag
116 d = n*(n + 1)//2
118 # Dual variable.
119 vars.append((SymmetricVariable.make_var_type(dim=d, bnd=0), k))
121 # Conic membership constraint.
122 cons.append((LMIConstraint.make_type(diag=n), k))
123 else:
124 assert False, "Unexpected constraint type."
126 # HACK: See above.
127 vars.append((RealVariable.make_var_type(dim=1, bnd=2), 1))
129 # Option changes.
130 nd_opts = footprint.options.updated(
131 dualize=False,
132 primals=footprint.options.duals,
133 duals=footprint.options.primals
134 ).nondefaults
136 return Footprint.from_types(dir, obj, vars, cons, nd_opts)
138 def forward(self):
139 """Implement :meth:`~.reformulation.Reformulation.forward`."""
140 self._pc2dv = {} # Maps primal constraints to dual variables.
141 self._pv2dc = {} # Maps primal variables to dual constraints.
143 P = self.input
145 # Create the dual problem from scratch.
146 D = self.output = Problem(copyOptions=P.options)
148 # Reset the 'dualize' option so that solvers can load the dual.
149 D.options.dualize = False
151 # Swap 'primals' and 'duals' option.
152 D.options.primals = P.options.duals
153 D.options.duals = P.options.primals
155 # Retrieve the primal objective in standard form.
156 original_primal_direction, primal_objective = P.objective.normalized
157 if original_primal_direction == "max":
158 primal_objective = -primal_objective
160 # Start building the dual objective function.
161 dual_objective = primal_objective.cst
163 # Start building the dual linear equalities for each primal variable.
164 obj_coefs = primal_objective._linear_coefs
165 dual_equality_terms = {var: -Constant(obj_coefs[var].T)
166 if var in obj_coefs else AffineExpression.zero(var.dim)
167 for var in P.variables.values()}
169 # Turn variable bounds into affine inequality constraints.
170 # NOTE: If bound-free variables are required by another reformulation or
171 # solver in the future, there should be a reformulation for this.
172 bound_cons = []
173 for primal_var in P.variables.values():
174 bound_constraint = primal_var.bound_constraint
175 if bound_constraint:
176 bound_cons.append(bound_constraint)
178 # Handle the primal constraints.
179 # NOTE: Constraint conversion works as follows:
180 # 1. Transform constraint into conic standard form:
181 # fx = Ax-b >_K 0.
182 # 2. Create a dual variable for the constraint.
183 # 3. Constrain the dual variable to be a member of the dual
184 # cone K*.
185 # 4. Extend the dual's linear equalities for each primal variable.
186 # 5. Extend the dual's objective function.
187 # TODO: Consider adding a ConicConstraint abstract subclass of
188 # Constraint with a method conic_form that returns a pair of cone
189 # member (affine expression) and cone. The conic standard form
190 # function f(x) and the dual cone can then be obtained by negation
191 # and through a new Cone.dual method, respectively.
192 for primal_con in chain(P.constraints.values(), bound_cons):
193 if isinstance(primal_con, AffineConstraint):
194 # Transform constraint into conic standard form.
195 fx = primal_con.ge0.vec
197 # Create a dual variable for the constraint.
198 self._pc2dv[primal_con] = dual_var = RealVariable(
199 "aff_{}".format(primal_con.id), len(fx))
201 # Constrain the dual variable to be a member of the dual cone.
202 # NOTE: In the equality case, the dual cone is the real field.
203 if primal_con.is_inequality():
204 D.add_constraint(dual_var >= 0)
206 # Extend the dual's linear equalities for each primal variable.
207 for var, coef in fx._linear_coefs.items():
208 dual_equality_terms[var] += coef.T*dual_var
210 # Extend the dual's objective function.
211 dual_objective -= fx._constant_coef.T*dual_var
213 elif isinstance(primal_con, SOCConstraint):
214 # Transform constraint into conic standard form.
215 fx = (primal_con.ub // primal_con.ne.vec)
217 # Create a dual variable for the constraint.
218 self._pc2dv[primal_con] = dual_var = RealVariable(
219 "soc_{}".format(primal_con.id), len(fx))
221 # Constrain the dual variable to be a member of the dual cone.
222 D.add_constraint(dual_var << soc())
224 # Extend the dual's linear equalities for each primal variable.
225 for var, coef in fx._linear_coefs.items():
226 dual_equality_terms[var] += coef.T*dual_var
228 # Extend the dual's objective function.
229 dual_objective -= fx._constant_coef.T*dual_var
230 elif isinstance(primal_con, RSOCConstraint):
231 # Transform constraint into conic standard form.
232 fx = (primal_con.ub1 // primal_con.ub2 // primal_con.ne.vec)
234 # Create a dual variable for the constraint.
235 self._pc2dv[primal_con] = dual_var = RealVariable(
236 "rso_{}".format(primal_con.id), len(fx))
238 # Constrain the dual variable to be a member of the dual cone.
239 D.add_constraint(dual_var << rsoc(4))
241 # Extend the dual's linear equalities for each primal variable.
242 for var, coef in fx._linear_coefs.items():
243 dual_equality_terms[var] += coef.T*dual_var
245 # Extend the dual's objective function.
246 dual_objective -= fx._constant_coef.T*dual_var
247 elif isinstance(primal_con, LMIConstraint):
248 # Transform constraint into conic standard form.
249 fx = primal_con.psd
251 # Create a dual variable for the constraint.
252 self._pc2dv[primal_con] = dual_var = SymmetricVariable(
253 "lmi_{}".format(primal_con.id), fx.shape)
255 # Constrain the dual variable to be a member of the dual cone.
256 D.add_constraint(dual_var >> 0)
258 # Extend the dual's linear equalities for each primal variable.
259 for var, coef in fx._linear_coefs.items():
260 dual_equality_terms[var] += coef.T*dual_var.vec
262 # Extend the dual's objective function.
263 dual_objective -= fx._constant_coef.T*dual_var.vec
264 else:
265 assert False, "Unexpected constraint type."
267 # Add the finished linear equalities for each variable.
268 for var, term in dual_equality_terms.items():
269 self._pv2dc[var] = D.add_constraint(term == 0)
271 # Adjust dual optimization direction to obtain the duality property.
272 if original_primal_direction == "max":
273 dual_direction = "min"
274 dual_objective = -dual_objective
275 else:
276 dual_direction = "max"
278 # HACK: Add a zero-fixed variable to the objective so that it is never
279 # constant. This makes prediction succeed in any case.
280 dual_objective += RealVariable("zero", 1, lower=0, upper=0)
282 # Set the finished dual objective.
283 D.objective = dual_direction, dual_objective
285 def update(self):
286 """Implement :meth:`~.reformulation.Reformulation.update`."""
287 # TODO: Implement updates to an existing dualization. This is optional.
288 raise NotImplementedError
290 def backward(self, solution):
291 """Implement :meth:`~.reformulation.Reformulation.backward`."""
292 P = self.input
294 # Require primal values to be properly vectorized.
295 if not solution.vectorizedPrimals:
296 raise NotImplementedError(
297 "Solving with dualize=True is not supported with solvers that "
298 "produce unvectorized primal solutions.")
300 # Retrieve primals for the primal problem.
301 primals = {}
302 for var in P.variables.values():
303 con = self._pv2dc[var]
305 if con in solution.duals and solution.duals[con] is not None:
306 primal = -solution.duals[con]
307 else:
308 primal = None
310 if primal is not None:
311 primals[var] = cvxopt.matrix(primal)
312 else:
313 primals[var] = None
315 # Retrieve duals for the primal problem.
316 duals = {}
317 for con in P.constraints.values():
318 var = self._pc2dv[con]
319 dual = solution.primals[var] if var in solution.primals else None
321 if dual is not None:
322 duals[con] = var._vec.devectorize(cvxopt.matrix(dual))
323 else:
324 duals[con] = None
326 # Swap infeasible/unbounded for a problem status.
327 if solution.problemStatus == PS_UNBOUNDED:
328 problemStatus = PS_INFEASIBLE
329 elif solution.problemStatus == PS_INFEASIBLE:
330 problemStatus = PS_UNBOUNDED
331 else:
332 problemStatus = solution.problemStatus
334 return Solution(
335 primals=primals,
336 duals=duals,
337 problem=solution.problem,
338 solver=solution.solver,
339 primalStatus=solution.dualStatus, # Intentional swap.
340 dualStatus=solution.primalStatus, # Intentional swap.
341 problemStatus=problemStatus,
342 searchTime=solution.searchTime,
343 info=solution.info,
344 vectorizedPrimals=True)
347# --------------------------------------
348__all__ = api_end(_API_START, globals())