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# ------------------------------------------------------------------------------ 

18 

19"""Implementation of :class:`Dualization`.""" 

20 

21from itertools import chain 

22 

23import cvxopt 

24 

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 

36 

37_API_START = api_start(globals()) 

38# ------------------------------- 

39 

40 

41class Dualization(Reformulation): 

42 """Lagrange dual problem reformulation.""" 

43 

44 SUPPORTED = Specification( 

45 objectives=[ 

46 AffineExpression], 

47 variables=CONTINUOUS_VARTYPES, 

48 constraints=[ 

49 DummyConstraint, 

50 AffineConstraint, 

51 SOCConstraint, 

52 RSOCConstraint, 

53 LMIConstraint]) 

54 

55 @classmethod 

56 def supports(cls, footprint): 

57 """Implement :meth:`~.reformulation.Reformulation.supports`.""" 

58 return footprint.options.dualize and footprint in cls.SUPPORTED 

59 

60 @classmethod 

61 def predict(cls, footprint): 

62 """Implement :meth:`~.reformulation.Reformulation.predict`.""" 

63 vars = [] 

64 cons = [] 

65 

66 # Objective direction. 

67 dir = "min" if footprint.direction == "max" else "max" 

68 

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) 

74 

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)) 

80 

81 if vartype.subtype.bnd: 

82 bound_cons.append((AffineConstraint.make_type( 

83 dim=vartype.subtype.bnd, eq=False), k)) 

84 

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 

91 

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 

98 

99 # Dual variable. 

100 vars.append((RealVariable.make_var_type(dim=d, bnd=0), k)) 

101 

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 

107 

108 # Dual variable. 

109 vars.append((RealVariable.make_var_type(dim=d, bnd=0), k)) 

110 

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 

116 

117 # Dual variable. 

118 vars.append((SymmetricVariable.make_var_type(dim=d, bnd=0), k)) 

119 

120 # Conic membership constraint. 

121 cons.append((LMIConstraint.make_type(diag=n), k)) 

122 else: 

123 assert False, "Unexpected constraint type." 

124 

125 # HACK: See above. 

126 vars.append((RealVariable.make_var_type(dim=1, bnd=2), 1)) 

127 

128 # Option changes. 

129 nd_opts = footprint.options.updated( 

130 dualize=False, 

131 primals=footprint.options.duals, 

132 duals=footprint.options.primals 

133 ).nondefaults 

134 

135 return Footprint.from_types(dir, obj, vars, cons, nd_opts) 

136 

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. 

141 

142 P = self.input 

143 

144 # Create the dual problem from scratch. 

145 D = self.output = Problem(copyOptions=P.options) 

146 

147 # Reset the 'dualize' option so that solvers can load the dual. 

148 D.options.dualize = False 

149 

150 # Swap 'primals' and 'duals' option. 

151 D.options.primals = P.options.duals 

152 D.options.duals = P.options.primals 

153 

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 

158 

159 # Start building the dual objective function. 

160 dual_objective = primal_objective.cst 

161 

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()} 

167 

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) 

176 

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 

197 

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) 

203 

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 

207 

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) 

213 

214 # Create a dual variable for the constraint. 

215 self._pc2dv[primal_con] = dual_var = RealVariable( 

216 "soc_{}".format(primal_con.id), len(fx)) 

217 

218 # Constrain the dual variable to be a member of the dual cone. 

219 D.add_constraint(dual_var << soc()) 

220 

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 

224 

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) 

230 

231 # Create a dual variable for the constraint. 

232 self._pc2dv[primal_con] = dual_var = RealVariable( 

233 "rso_{}".format(primal_con.id), len(fx)) 

234 

235 # Constrain the dual variable to be a member of the dual cone. 

236 D.add_constraint(dual_var << rsoc(4)) 

237 

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 

241 

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 

247 

248 # Create a dual variable for the constraint. 

249 self._pc2dv[primal_con] = dual_var = SymmetricVariable( 

250 "lmi_{}".format(primal_con.id), fx.shape) 

251 

252 # Constrain the dual variable to be a member of the dual cone. 

253 D.add_constraint(dual_var >> 0) 

254 

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 

258 

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." 

263 

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) 

267 

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" 

274 

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) 

278 

279 # Set the finished dual objective. 

280 D.objective = dual_direction, dual_objective 

281 

282 def update(self): 

283 """Implement :meth:`~.reformulation.Reformulation.update`.""" 

284 # TODO: Implement updates to an existing dualization. This is optional. 

285 raise NotImplementedError 

286 

287 def backward(self, solution): 

288 """Implement :meth:`~.reformulation.Reformulation.backward`.""" 

289 P = self.input 

290 

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.") 

296 

297 # Retrieve primals for the primal problem. 

298 primals = {} 

299 for var in P.variables.values(): 

300 con = self._pv2dc[var] 

301 

302 if con in solution.duals and solution.duals[con] is not None: 

303 primal = -solution.duals[con] 

304 else: 

305 primal = None 

306 

307 if primal is not None: 

308 primals[var] = cvxopt.matrix(primal) 

309 else: 

310 primals[var] = None 

311 

312 # Retrieve duals for the primal problem. 

313 duals = {} 

314 for con in P.constraints.values(): 

315 if isinstance(con, DummyConstraint): 

316 continue 

317 

318 var = self._pc2dv[con] 

319 dual = solution.primals[var] if var in solution.primals else None 

320 

321 if dual is not None: 

322 duals[con] = var._vec.devectorize(cvxopt.matrix(dual)) 

323 else: 

324 duals[con] = None 

325 

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 

333 

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) 

346 

347 

348# -------------------------------------- 

349__all__ = api_end(_API_START, globals())