Coverage for picos/reforms/reform_dualize.py: 97.90%
143 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
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 with conic membership as bounds.
93 vars.append((RealVariable.make_var_type(
94 dim=d, bnd=0 if contype.subtype.eq else d), k))
95 elif contype.clstype is SOCConstraint:
96 n = contype.subtype.argdim
97 d = n + 1
99 # Dual variable.
100 vars.append((RealVariable.make_var_type(dim=d, bnd=0), k))
102 # Conic membership constraint.
103 cons.append((SOCConstraint.make_type(argdim=n), k))
104 elif contype.clstype is RSOCConstraint:
105 n = contype.subtype.argdim
106 d = n + 2
108 # Dual variable.
109 vars.append((RealVariable.make_var_type(dim=d, bnd=0), k))
111 # Conic membership constraint.
112 cons.append((RSOCConstraint.make_type(argdim=n), k))
113 elif contype.clstype is LMIConstraint:
114 n = contype.subtype.diag
115 d = n*(n + 1)//2
117 # Dual variable.
118 vars.append((SymmetricVariable.make_var_type(dim=d, bnd=0), k))
120 # Conic membership constraint.
121 cons.append((LMIConstraint.make_type(diag=n), k))
122 else:
123 assert False, "Unexpected constraint type."
125 # HACK: See above.
126 vars.append((RealVariable.make_var_type(dim=1, bnd=2), 1))
128 # Option changes.
129 nd_opts = footprint.options.updated(
130 dualize=False,
131 primals=footprint.options.duals,
132 duals=footprint.options.primals
133 ).nondefaults
135 return Footprint.from_types(dir, obj, vars, cons, nd_opts)
137 def forward(self):
138 """Implement :meth:`~.reformulation.Reformulation.forward`."""
139 self._pc2dv = {} # Maps primal constraints to dual variables.
140 self._pv2dc = {} # Maps primal variables to dual constraints.
142 P = self.input
144 # Create the dual problem from scratch.
145 D = self.output = Problem(copyOptions=P.options)
147 # Reset the 'dualize' option so that solvers can load the dual.
148 D.options.dualize = False
150 # Swap 'primals' and 'duals' option.
151 D.options.primals = P.options.duals
152 D.options.duals = P.options.primals
154 # Retrieve the primal objective in standard form.
155 original_primal_direction, primal_objective = P.objective.normalized
156 if original_primal_direction == "max":
157 primal_objective = -primal_objective
159 # Start building the dual objective function.
160 dual_objective = primal_objective.cst
162 # Start building the dual linear equalities for each primal variable.
163 obj_coefs = primal_objective._linear_coefs
164 dual_equality_terms = {var: -Constant(obj_coefs[var].T)
165 if var in obj_coefs else AffineExpression.zero(var.dim)
166 for var in P.variables.values()}
168 # Turn variable bounds into affine inequality constraints.
169 # NOTE: If bound-free variables are required by another reformulation or
170 # solver in the future, there should be a reformulation for this.
171 bound_cons = []
172 for primal_var in P.variables.values():
173 bound_constraint = primal_var.bound_constraint
174 if bound_constraint:
175 bound_cons.append(bound_constraint)
177 # Handle the primal constraints.
178 # NOTE: Constraint conversion works as follows:
179 # 1. Transform constraint into conic standard form:
180 # fx = Ax-b >_K 0.
181 # 2. Create a dual variable for the constraint.
182 # 3. Constrain the dual variable to be a member of the dual
183 # cone K*.
184 # 4. Extend the dual's linear equalities for each primal variable.
185 # 5. Extend the dual's objective function.
186 # TODO: Consider adding a ConicConstraint abstract subclass of
187 # Constraint with a method conic_form that returns a pair of cone
188 # member (affine expression) and cone. The conic standard form
189 # function f(x) and the dual cone can then be obtained by negation
190 # and through a new Cone.dual method, respectively.
191 for primal_con in chain(P.constraints.values(), bound_cons):
192 if isinstance(primal_con, DummyConstraint):
193 pass
194 elif isinstance(primal_con, AffineConstraint):
195 # Transform constraint into conic standard form.
196 fx = primal_con.ge0.vec
198 # Create a dual variable for the constraint.
199 # NOTE: Conic membership is expressed through variable bounds.
200 self._pc2dv[primal_con] = dual_var = RealVariable(
201 "aff_{}".format(primal_con.id), len(fx),
202 lower=0 if primal_con.is_inequality() else None)
204 # Extend the dual's linear equalities for each primal variable.
205 for var, coef in fx._linear_coefs.items():
206 dual_equality_terms[var] += coef.T*dual_var
208 # Extend the dual's objective function.
209 dual_objective -= fx._constant_coef.T*dual_var
210 elif isinstance(primal_con, SOCConstraint):
211 # Transform constraint into conic standard form.
212 fx = (primal_con.ub // primal_con.ne.vec)
214 # Create a dual variable for the constraint.
215 self._pc2dv[primal_con] = dual_var = RealVariable(
216 "soc_{}".format(primal_con.id), len(fx))
218 # Constrain the dual variable to be a member of the dual cone.
219 D.add_constraint(dual_var << soc())
221 # Extend the dual's linear equalities for each primal variable.
222 for var, coef in fx._linear_coefs.items():
223 dual_equality_terms[var] += coef.T*dual_var
225 # Extend the dual's objective function.
226 dual_objective -= fx._constant_coef.T*dual_var
227 elif isinstance(primal_con, RSOCConstraint):
228 # Transform constraint into conic standard form.
229 fx = (primal_con.ub1 // primal_con.ub2 // primal_con.ne.vec)
231 # Create a dual variable for the constraint.
232 self._pc2dv[primal_con] = dual_var = RealVariable(
233 "rso_{}".format(primal_con.id), len(fx))
235 # Constrain the dual variable to be a member of the dual cone.
236 D.add_constraint(dual_var << rsoc(4))
238 # Extend the dual's linear equalities for each primal variable.
239 for var, coef in fx._linear_coefs.items():
240 dual_equality_terms[var] += coef.T*dual_var
242 # Extend the dual's objective function.
243 dual_objective -= fx._constant_coef.T*dual_var
244 elif isinstance(primal_con, LMIConstraint):
245 # Transform constraint into conic standard form.
246 fx = primal_con.psd
248 # Create a dual variable for the constraint.
249 self._pc2dv[primal_con] = dual_var = SymmetricVariable(
250 "lmi_{}".format(primal_con.id), fx.shape)
252 # Constrain the dual variable to be a member of the dual cone.
253 D.add_constraint(dual_var >> 0)
255 # Extend the dual's linear equalities for each primal variable.
256 for var, coef in fx._linear_coefs.items():
257 dual_equality_terms[var] += coef.T*dual_var.vec
259 # Extend the dual's objective function.
260 dual_objective -= fx._constant_coef.T*dual_var.vec
261 else:
262 assert False, "Unexpected constraint type."
264 # Add the finished linear equalities for each variable.
265 for var, term in dual_equality_terms.items():
266 self._pv2dc[var] = D.add_constraint(term == 0)
268 # Adjust dual optimization direction to obtain the duality property.
269 if original_primal_direction == "max":
270 dual_direction = "min"
271 dual_objective = -dual_objective
272 else:
273 dual_direction = "max"
275 # HACK: Add a zero-fixed variable to the objective so that it is never
276 # constant. This makes prediction succeed in any case.
277 dual_objective += RealVariable("zero", 1, lower=0, upper=0)
279 # Set the finished dual objective.
280 D.objective = dual_direction, dual_objective
282 def update(self):
283 """Implement :meth:`~.reformulation.Reformulation.update`."""
284 # TODO: Implement updates to an existing dualization. This is optional.
285 raise NotImplementedError
287 def backward(self, solution):
288 """Implement :meth:`~.reformulation.Reformulation.backward`."""
289 P = self.input
291 # Require primal values to be properly vectorized.
292 if not solution.vectorizedPrimals:
293 raise NotImplementedError(
294 "Solving with dualize=True is not supported with solvers that "
295 "produce unvectorized primal solutions.")
297 # Retrieve primals for the primal problem.
298 primals = {}
299 for var in P.variables.values():
300 con = self._pv2dc[var]
302 if con in solution.duals and solution.duals[con] is not None:
303 primal = -solution.duals[con]
304 else:
305 primal = None
307 if primal is not None:
308 primals[var] = cvxopt.matrix(primal)
309 else:
310 primals[var] = None
312 # Retrieve duals for the primal problem.
313 duals = {}
314 for con in P.constraints.values():
315 if isinstance(con, DummyConstraint):
316 continue
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,
345 reportedValue=solution.reportedValue)
348# --------------------------------------
349__all__ = api_end(_API_START, globals())