Hide keyboard shortcuts

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

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. 

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

94 

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 

102 

103 # Dual variable. 

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

105 

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 

111 

112 # Dual variable. 

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

114 

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 

120 

121 # Dual variable. 

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

123 

124 # Conic membership constraint. 

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

126 else: 

127 assert False, "Unexpected constraint type." 

128 

129 # HACK: See above. 

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

131 

132 # Option changes. 

133 nd_opts = footprint.options.updated( 

134 dualize=False, 

135 primals=footprint.options.duals, 

136 duals=footprint.options.primals 

137 ).nondefaults 

138 

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

140 

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. 

145 

146 P = self.input 

147 

148 # Create the dual problem from scratch. 

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

150 

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

152 D.options.dualize = False 

153 

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

155 D.options.primals = P.options.duals 

156 D.options.duals = P.options.primals 

157 

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 

162 

163 # Start building the dual objective function. 

164 dual_objective = primal_objective.cst 

165 

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

171 

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) 

180 

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 

201 

202 # Create a dual variable for the constraint. 

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

204 "aff_{}".format(primal_con.id), len(fx)) 

205 

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) 

210 

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 

214 

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) 

220 

221 # Create a dual variable for the constraint. 

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

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

224 

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

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

227 

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 

231 

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) 

237 

238 # Create a dual variable for the constraint. 

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

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

241 

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

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

244 

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 

248 

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 

254 

255 # Create a dual variable for the constraint. 

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

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

258 

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

260 D.add_constraint(dual_var >> 0) 

261 

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 

265 

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

270 

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) 

274 

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" 

281 

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) 

285 

286 # Set the finished dual objective. 

287 D.objective = dual_direction, dual_objective 

288 

289 def update(self): 

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

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

292 raise NotImplementedError 

293 

294 def backward(self, solution): 

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

296 P = self.input 

297 

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

303 

304 # Retrieve primals for the primal problem. 

305 primals = {} 

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

307 con = self._pv2dc[var] 

308 

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

310 primal = -solution.duals[con] 

311 else: 

312 primal = None 

313 

314 if primal is not None: 

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

316 else: 

317 primals[var] = None 

318 

319 # Retrieve duals for the primal problem. 

320 duals = {} 

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

322 if isinstance(con, DummyConstraint): 

323 continue 

324 

325 var = self._pc2dv[con] 

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

327 

328 if dual is not None: 

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

330 else: 

331 duals[con] = None 

332 

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 

340 

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) 

353 

354 

355# -------------------------------------- 

356__all__ = api_end(_API_START, globals())