Coverage for gpkit/constraints/sgp.py: 81%

Shortcuts on this page

r m x   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

169 statements  

1"""Implement the SequentialGeometricProgram class""" 

2import warnings as pywarnings 

3from time import time 

4from collections import defaultdict 

5import numpy as np 

6from ..exceptions import (InvalidGPConstraint, Infeasible, UnnecessarySGP, 

7 InvalidPosynomial, InvalidSGPConstraint) 

8from ..keydict import KeyDict 

9from ..nomials import Variable 

10from .gp import GeometricProgram 

11from ..nomials import PosynomialInequality, Posynomial 

12from .. import NamedVariables 

13from ..small_scripts import appendsolwarning, initsolwarning 

14 

15 

16EPS = 1e-6 # 1 +/- this is used in a few relative differences 

17 

18# pylint: disable=too-many-instance-attributes 

19class SequentialGeometricProgram: 

20 """Prepares a collection of signomials for a SP solve. 

21 

22 Arguments 

23 --------- 

24 cost : Posynomial 

25 Objective to minimize when solving 

26 constraints : list of Constraint or SignomialConstraint objects 

27 Constraints to maintain when solving (implicitly Signomials <= 1) 

28 verbosity : int (optional) 

29 Currently has no effect: SequentialGeometricPrograms don't know 

30 anything new after being created, unlike GeometricPrograms. 

31 

32 Attributes with side effects 

33 ---------------------------- 

34 `gps` is set during a solve 

35 `result` is set at the end of a solve 

36 

37 Examples 

38 -------- 

39 >>> gp = gpkit.geometric_program.SequentialGeometricProgram( 

40 # minimize 

41 x, 

42 [ # subject to 

43 1/x - y/x, # <= 1, implicitly 

44 y/10 # <= 1 

45 ]) 

46 >>> gp.solve() 

47 """ 

48 gps = solver_outs = _results = result = model = None 

49 with NamedVariables("RelaxPCCP"): 

50 slack = Variable("C") 

51 

52 def __init__(self, cost, model, substitutions, 

53 *, use_pccp=True, pccp_penalty=2e2, **kwargs): 

54 self.cost = cost 

55 self.pccp_penalty = pccp_penalty 

56 if cost.any_nonpositive_cs: 

57 raise InvalidPosynomial("""an SGP's cost must be Posynomial 

58 

59 The equivalent of a Signomial objective can be constructed by constraining 

60 a dummy variable `z` to be greater than the desired Signomial objective `s` 

61 (z >= s) and then minimizing that dummy variable.""") 

62 self.gpconstraints, self.sgpconstraints = [], [] 

63 if not use_pccp: 

64 self.slack = 1 

65 else: 

66 self.gpconstraints.append(self.slack >= 1) 

67 cost *= self.slack**pccp_penalty 

68 self.approxconstraints = [] 

69 self.sgpvks = set() 

70 x0 = KeyDict(substitutions) 

71 x0.vks = model.vks # for string access and so forth 

72 for cs in model.flat(): 

73 try: 

74 if not hasattr(cs, "as_hmapslt1"): 

75 raise InvalidGPConstraint(cs) 

76 if not isinstance(cs, PosynomialInequality): 

77 cs.as_hmapslt1(substitutions) # gp-compatible? 

78 self.gpconstraints.append(cs) 

79 except InvalidGPConstraint: 

80 if not hasattr(cs, "as_gpconstr"): 

81 raise InvalidSGPConstraint(cs) 

82 self.sgpconstraints.append(cs) 

83 for hmaplt1 in cs.as_gpconstr(x0).as_hmapslt1({}): 

84 constraint = (Posynomial(hmaplt1) <= self.slack) 

85 constraint.generated_by = cs 

86 self.approxconstraints.append(constraint) 

87 self.sgpvks.update(constraint.vks) 

88 if not self.sgpconstraints: 

89 raise UnnecessarySGP("""Model valid as a Geometric Program. 

90 

91SequentialGeometricPrograms should only be created with Models containing 

92Signomial Constraints, since Models without Signomials have global 

93solutions and can be solved with 'Model.solve()'.""") 

94 self._gp = GeometricProgram( 

95 cost, self.approxconstraints + self.gpconstraints, 

96 substitutions, **kwargs) 

97 self._gp.x0 = x0 

98 self.a_idxs = defaultdict(list) 

99 last_cost_mon = self._gp.k[0] 

100 first_gp_mon = sum(self._gp.k[:1+len(self.approxconstraints)]) 

101 for row_idx, m_idx in enumerate(self._gp.A.row): 

102 if last_cost_mon <= m_idx <= first_gp_mon: 

103 self.a_idxs[self._gp.p_idxs[m_idx]].append(row_idx) 

104 

105 # pylint: disable=too-many-locals,too-many-branches,too-many-statements 

106 def localsolve(self, solver=None, *, verbosity=1, x0=None, reltol=1e-4, 

107 iteration_limit=50, err_on_relax=True, **solveargs): 

108 """Locally solves a SequentialGeometricProgram and returns the solution. 

109 

110 Arguments 

111 --------- 

112 solver : str or function (optional) 

113 By default uses one of the solvers found during installation. 

114 If set to "mosek", "mosek_cli", or "cvxopt", uses that solver. 

115 If set to a function, passes that function cs, A, p_idxs, and k. 

116 verbosity : int (optional) 

117 If greater than 0, prints solve time and number of iterations. 

118 Each GP is created and solved with verbosity one less than this, so 

119 if greater than 1, prints solver name and time for each GP. 

120 x0 : dict (optional) 

121 Initial location to approximate signomials about. 

122 reltol : float 

123 Iteration ends when this is greater than the distance between two 

124 consecutive solve's objective values. 

125 iteration_limit : int 

126 Maximum GP iterations allowed. 

127 mutategp: boolean 

128 Prescribes whether to mutate the previously generated GP 

129 or to create a new GP with every solve. 

130 **solveargs : 

131 Passed to solver function. 

132 

133 Returns 

134 ------- 

135 result : dict 

136 A dictionary containing the translated solver result. 

137 """ 

138 self.gps, self.solver_outs, self._results = [], [], [] 

139 starttime = time() 

140 if verbosity > 0: 

141 print("Starting a sequence of GP solves") 

142 print(" for %i free variables" % len(self.sgpvks)) 

143 print(" in %i locally-GP constraints" % len(self.sgpconstraints)) 

144 print(" and for %i free variables" % len(self._gp.varlocs)) 

145 print(" in %i posynomial inequalities." % len(self._gp.k)) 

146 prevcost, cost, rel_improvement = None, None, None 

147 while rel_improvement is None or rel_improvement > reltol: 

148 prevcost = cost 

149 if len(self.gps) > iteration_limit: 

150 raise Infeasible( 

151 "Unsolved after %s iterations. Check `m.program.results`;" 

152 " if they're converging, try `.localsolve(...," 

153 " iteration_limit=NEWLIMIT)`." % len(self.gps)) 

154 gp = self.gp(x0, cleanx0=(len(self.gps) >= 1)) # clean the first x0 

155 self.gps.append(gp) # NOTE: SIDE EFFECTS 

156 if verbosity > 1: 

157 print("\nGP Solve %i" % len(self.gps)) 

158 if verbosity > 2: 

159 print("===============") 

160 solver_out = gp.solve(solver, verbosity=verbosity-1, 

161 gen_result=False, **solveargs) 

162 self.solver_outs.append(solver_out) 

163 cost = float(solver_out["objective"]) 

164 x0 = dict(zip(gp.varlocs, np.exp(solver_out["primal"]))) 

165 if verbosity > 2: 

166 result = gp.generate_result(solver_out, verbosity=verbosity-3) 

167 self._results.append(result) 

168 vartable = result.table(self.sgpvks, tables=["freevariables"]) 

169 vartable = "\n" + vartable.replace("Free", "SGP", 1) 

170 print(vartable) 

171 elif verbosity > 1: 

172 print("Solved cost was %.4g." % cost) 

173 if prevcost is None: 

174 continue 

175 rel_improvement = (prevcost - cost)/(prevcost + cost) 

176 if cost/prevcost >= 1 + 10*EPS: 

177 pywarnings.warn( 

178 "SGP not convergent: Cost rose by %.2g%% (%.6g to %.6g) on" 

179 " GP solve %i. Details can be found in `m.program.results`" 

180 " or by solving at a higher verbosity. Note convergence" 

181 " is not guaranteed for models with SignomialEqualities." 

182 % (100*(cost - prevcost)/prevcost, 

183 prevcost, cost, len(self.gps))) 

184 rel_improvement = cost = None 

185 # solved successfully! 

186 self.result = gp.generate_result(solver_out, verbosity=verbosity-3) 

187 self.result["soltime"] = time() - starttime 

188 if verbosity > 1: 

189 print() 

190 if verbosity > 0: 

191 print("Solving took %.3g seconds and %i GP solves." 

192 % (self.result["soltime"], len(self.gps))) 

193 if hasattr(self.slack, "key"): 

194 initsolwarning(self.result, "Slack Non-GP Constraints") 

195 excess_slack = self.result["variables"][self.slack.key] - 1 # pylint: disable=no-member 

196 if excess_slack > EPS: 

197 msg = ("Final PCCP solution let non-GP constraints slacken by" 

198 " %.2g%%." % (100*excess_slack)) 

199 expl = (msg + 

200 " Calling .localsolve(pccp_penalty=...) with a higher" 

201 " `pccp_penalty` (it was %.3g this time) will reduce" 

202 " slack if the model is solvable with less. To verify" 

203 " that the slack is needed, generate an SGP with" 

204 " `use_pccp=False` and start it from this model's" 

205 " solution: e.g. `m.localsolve(use_pccp=False, x0=" 

206 "m.solution[\"variables\"])`." % self.pccp_penalty) 

207 if err_on_relax: 

208 raise Infeasible(expl) 

209 appendsolwarning(msg, (1 + excess_slack), self.result, 

210 "Slack Non-GP Constraints") 

211 if verbosity > -1: 

212 print(expl) 

213 self.result["cost function"] = self.cost 

214 del self.result["freevariables"][self.slack.key] # pylint: disable=no-member 

215 del self.result["variables"][self.slack.key] # pylint: disable=no-member 

216 if "sensitivities" in self.result: # not true for MIGP 

217 del self.result["sensitivities"]["variables"][self.slack.key] # pylint: disable=no-member 

218 del self.result["sensitivities"]["variablerisk"][self.slack.key] # pylint: disable=no-member 

219 slcon = self.gpconstraints[0] 

220 slconsenss = self.result["sensitivities"]["constraints"][slcon] 

221 del self.result["sensitivities"]["constraints"][slcon] 

222 # TODO: create constraint in RelaxPCCP namespace 

223 self.result["sensitivities"]["models"][""] -= slconsenss 

224 if not self.result["sensitivities"]["models"][""]: 

225 del self.result["sensitivities"]["models"][""] 

226 return self.result 

227 

228 @property 

229 def results(self): 

230 "Creates and caches results from the raw solver_outs" 

231 if not self._results: 

232 self._results = [gp.generate_result(s_o, dual_check=False) 

233 for gp, s_o in zip(self.gps, self.solver_outs)] 

234 return self._results 

235 

236 def gp(self, x0=None, *, cleanx0=False): 

237 "Update self._gp for x0 and return it." 

238 if not x0: 

239 return self._gp # return last generated 

240 if not cleanx0: 

241 cleanedx0 = KeyDict() 

242 cleanedx0.vks = self._gp.x0.vks 

243 cleanedx0.update(x0) 

244 x0 = cleanedx0 

245 self._gp.x0.update({vk: x0[vk] for vk in self.sgpvks if vk in x0}) 

246 p_idx = 0 

247 for sgpc in self.sgpconstraints: 

248 for hmaplt1 in sgpc.as_gpconstr(self._gp.x0).as_hmapslt1({}): 

249 approxc = self.approxconstraints[p_idx] 

250 approxc.left = self.slack 

251 approxc.right.hmap = hmaplt1 

252 approxc.unsubbed = [Posynomial(hmaplt1)/self.slack] 

253 p_idx += 1 # p_idx=0 is the cost; sp constraints are after it 

254 hmap, = approxc.as_hmapslt1(self._gp.substitutions) 

255 self._gp.hmaps[p_idx] = hmap 

256 m_idx = self._gp.m_idxs[p_idx].start 

257 a_idxs = list(self.a_idxs[p_idx]) # A's entries we can modify 

258 for i, (exp, c) in enumerate(hmap.items()): 

259 self._gp.exps[m_idx + i] = exp 

260 self._gp.cs[m_idx + i] = c 

261 for var, x in exp.items(): 

262 try: # modify a particular A entry 

263 row_idx = a_idxs.pop() 

264 self._gp.A.row[row_idx] = m_idx + i 

265 self._gp.A.col[row_idx] = self._gp.varidxs[var] 

266 self._gp.A.data[row_idx] = x 

267 except IndexError: # numbers of exps increased 

268 self.a_idxs[p_idx].append(len(self._gp.A.row)) 

269 self._gp.A.row.append(m_idx + i) 

270 self._gp.A.col.append(self._gp.varidxs[var]) 

271 self._gp.A.data.append(x) 

272 for row_idx in a_idxs: # number of exps decreased 

273 self._gp.A.row[row_idx] = 0 # zero out this entry 

274 self._gp.A.col[row_idx] = 0 

275 self._gp.A.data[row_idx] = 0 

276 return self._gp