Coverage for gpkit/constraints/sgp.py: 82%
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
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
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
16EPS = 1e-6 # 1 +/- this is used in a few relative differences
18# pylint: disable=too-many-instance-attributes
19class SequentialGeometricProgram:
20 """Prepares a collection of signomials for a SP solve.
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.
32 Attributes with side effects
33 ----------------------------
34 `gps` is set during a solve
35 `result` is set at the end of a solve
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")
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
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.
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)
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, **solveargs):
108 """Locally solves a SequentialGeometricProgram and returns the solution.
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.
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 appendsolwarning(msg, (1 + excess_slack), self.result,
200 "Slack Non-GP Constraints")
201 if verbosity > -1:
202 print(msg +
203 " Calling .localsolve(pccp_penalty=...) with a higher"
204 " `pccp_penalty` (it was %.3g this time) will reduce"
205 " slack if the model is solvable with less. To verify"
206 " that the slack is needed, generate an SGP with"
207 " `use_pccp=False` and start it from this model's"
208 " solution: e.g. `m.localsolve(use_pccp=False, x0="
209 "m.solution[\"variables\"])`." % self.pccp_penalty)
210 self.result["cost function"] = self.cost
211 del self.result["freevariables"][self.slack.key] # pylint: disable=no-member
212 del self.result["variables"][self.slack.key] # pylint: disable=no-member
213 del self.result["sensitivities"]["variables"][self.slack.key] # pylint: disable=no-member
214 del self.result["sensitivities"]["variablerisk"][self.slack.key] # pylint: disable=no-member
215 slcon = self.gpconstraints[0]
216 slconsenss = self.result["sensitivities"]["constraints"][slcon]
217 del self.result["sensitivities"]["constraints"][slcon]
218 # TODO: create constraint in RelaxPCCP namespace
219 self.result["sensitivities"]["models"][""] -= slconsenss
220 if not self.result["sensitivities"]["models"][""]:
221 del self.result["sensitivities"]["models"][""]
222 return self.result
224 @property
225 def results(self):
226 "Creates and caches results from the raw solver_outs"
227 if not self._results:
228 self._results = [gp.generate_result(s_o, dual_check=False)
229 for gp, s_o in zip(self.gps, self.solver_outs)]
230 return self._results
232 def gp(self, x0=None, *, cleanx0=False):
233 "Update self._gp for x0 and return it."
234 if not x0:
235 return self._gp # return last generated
236 if not cleanx0:
237 cleanedx0 = KeyDict()
238 cleanedx0.vks = self._gp.x0.vks
239 cleanedx0.update(x0)
240 x0 = cleanedx0
241 self._gp.x0.update({vk: x0[vk] for vk in self.sgpvks if vk in x0})
242 p_idx = 0
243 for sgpc in self.sgpconstraints:
244 for hmaplt1 in sgpc.as_gpconstr(self._gp.x0).as_hmapslt1({}):
245 approxc = self.approxconstraints[p_idx]
246 approxc.left = self.slack
247 approxc.right.hmap = hmaplt1
248 approxc.unsubbed = [Posynomial(hmaplt1)/self.slack]
249 p_idx += 1 # p_idx=0 is the cost; sp constraints are after it
250 hmap, = approxc.as_hmapslt1(self._gp.substitutions)
251 self._gp.hmaps[p_idx] = hmap
252 m_idx = self._gp.m_idxs[p_idx].start
253 a_idxs = list(self.a_idxs[p_idx]) # A's entries we can modify
254 for i, (exp, c) in enumerate(hmap.items()):
255 self._gp.exps[m_idx + i] = exp
256 self._gp.cs[m_idx + i] = c
257 for var, x in exp.items():
258 try: # modify a particular A entry
259 row_idx = a_idxs.pop()
260 self._gp.A.row[row_idx] = m_idx + i
261 self._gp.A.col[row_idx] = self._gp.varidxs[var]
262 self._gp.A.data[row_idx] = x
263 except IndexError: # numbers of exps increased
264 self.a_idxs[p_idx].append(len(self._gp.A.row))
265 self._gp.A.row.append(m_idx + i)
266 self._gp.A.col.append(self._gp.varidxs[var])
267 self._gp.A.data.append(x)
268 for row_idx in a_idxs: # number of exps decreased
269 self._gp.A.row[row_idx] = 0 # zero out this entry
270 self._gp.A.col[row_idx] = 0
271 self._gp.A.data[row_idx] = 0
272 return self._gp