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

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 AffineConstraint, 

50 SOCConstraint, 

51 RSOCConstraint, 

52 LMIConstraint]) 

53 

54 @classmethod 

55 def supports(cls, footprint): 

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

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

58 

59 @classmethod 

60 def predict(cls, footprint): 

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

62 vars = [] 

63 cons = [] 

64 

65 # Objective direction. 

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

67 

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) 

73 

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

79 

80 if vartype.subtype.bnd: 

81 bound_cons.append((AffineConstraint.make_type( 

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

83 

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 

88 

89 # Dual variable. 

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

91 

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 

99 

100 # Dual variable. 

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

102 

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 

108 

109 # Dual variable. 

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

111 

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 

117 

118 # Dual variable. 

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

120 

121 # Conic membership constraint. 

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

123 else: 

124 assert False, "Unexpected constraint type." 

125 

126 # HACK: See above. 

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

128 

129 # Option changes. 

130 nd_opts = footprint.options.updated( 

131 dualize=False, 

132 primals=footprint.options.duals, 

133 duals=footprint.options.primals 

134 ).nondefaults 

135 

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

137 

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. 

142 

143 P = self.input 

144 

145 # Create the dual problem from scratch. 

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

147 

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

149 D.options.dualize = False 

150 

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

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

153 D.options.duals = P.options.primals 

154 

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 

159 

160 # Start building the dual objective function. 

161 dual_objective = primal_objective.cst 

162 

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

168 

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) 

177 

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 

196 

197 # Create a dual variable for the constraint. 

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

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

200 

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) 

205 

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 

209 

210 # Extend the dual's objective function. 

211 dual_objective -= fx._constant_coef.T*dual_var 

212 

213 elif isinstance(primal_con, SOCConstraint): 

214 # Transform constraint into conic standard form. 

215 fx = (primal_con.ub // primal_con.ne.vec) 

216 

217 # Create a dual variable for the constraint. 

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

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

220 

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

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

223 

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 

227 

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) 

233 

234 # Create a dual variable for the constraint. 

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

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

237 

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

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

240 

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 

244 

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 

250 

251 # Create a dual variable for the constraint. 

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

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

254 

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

256 D.add_constraint(dual_var >> 0) 

257 

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 

261 

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

266 

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) 

270 

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" 

277 

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) 

281 

282 # Set the finished dual objective. 

283 D.objective = dual_direction, dual_objective 

284 

285 def update(self): 

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

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

288 raise NotImplementedError 

289 

290 def backward(self, solution): 

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

292 P = self.input 

293 

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

299 

300 # Retrieve primals for the primal problem. 

301 primals = {} 

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

303 con = self._pv2dc[var] 

304 

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

306 primal = -solution.duals[con] 

307 else: 

308 primal = None 

309 

310 if primal is not None: 

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

312 else: 

313 primals[var] = None 

314 

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 

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 

346 

347# -------------------------------------- 

348__all__ = api_end(_API_START, globals())