Coverage for gpkit/constraints/sgp.py: 80%
169 statements
« prev ^ index » next coverage.py v7.4.0, created at 2024-01-07 22:15 -0500
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 #pylint: disable=too-many-arguments, too-many-locals
53 def __init__(self, cost, model, substitutions,
54 *, use_pccp=True, pccp_penalty=2e2, **kwargs):
55 self.cost = cost
56 self.pccp_penalty = pccp_penalty
57 if cost.any_nonpositive_cs:
58 raise InvalidPosynomial("""an SGP's cost must be Posynomial
60 The equivalent of a Signomial objective can be constructed by constraining
61 a dummy variable `z` to be greater than the desired Signomial objective `s`
62 (z >= s) and then minimizing that dummy variable.""")
63 self.gpconstraints, self.sgpconstraints = [], []
64 if not use_pccp:
65 self.slack = 1
66 else:
67 self.gpconstraints.append(self.slack >= 1)
68 cost *= self.slack**pccp_penalty
69 self.approxconstraints = []
70 self.sgpvks = set()
71 x0 = KeyDict(substitutions)
72 x0.vks = model.vks # for string access and so forth
73 for cs in model.flat():
74 try:
75 if not hasattr(cs, "as_hmapslt1"):
76 raise InvalidGPConstraint(cs)
77 if not isinstance(cs, PosynomialInequality):
78 cs.as_hmapslt1(substitutions) # gp-compatible?
79 self.gpconstraints.append(cs)
80 except InvalidGPConstraint as err:
81 if not hasattr(cs, "as_gpconstr"):
82 raise InvalidSGPConstraint(cs) from err
83 self.sgpconstraints.append(cs)
84 for hmaplt1 in cs.as_gpconstr(x0).as_hmapslt1({}):
85 constraint = Posynomial(hmaplt1) <= self.slack
86 constraint.generated_by = cs
87 self.approxconstraints.append(constraint)
88 self.sgpvks.update(constraint.vks)
89 if not self.sgpconstraints:
90 raise UnnecessarySGP("""Model valid as a Geometric Program.
92SequentialGeometricPrograms should only be created with Models containing
93Signomial Constraints, since Models without Signomials have global
94solutions and can be solved with 'Model.solve()'.""")
95 self._gp = GeometricProgram(
96 cost, self.approxconstraints + self.gpconstraints,
97 substitutions, **kwargs)
98 self._gp.x0 = x0
99 self.a_idxs = defaultdict(list)
100 last_cost_mon = self._gp.k[0]
101 first_gp_mon = sum(self._gp.k[:1+len(self.approxconstraints)])
102 for row_idx, m_idx in enumerate(self._gp.A.row):
103 if last_cost_mon <= m_idx <= first_gp_mon:
104 self.a_idxs[self._gp.p_idxs[m_idx]].append(row_idx)
106 # pylint: disable=too-many-locals,too-many-branches,too-many-statements
107 # pylint: disable=too-many-arguments
108 def localsolve(self, solver=None, *, verbosity=1, x0=None, reltol=1e-4,
109 iteration_limit=50, err_on_relax=True, **solveargs):
110 """Locally solves a SequentialGeometricProgram and returns the solution.
112 Arguments
113 ---------
114 solver : str or function (optional)
115 By default uses one of the solvers found during installation.
116 If set to "mosek", "mosek_cli", or "cvxopt", uses that solver.
117 If set to a function, passes that function cs, A, p_idxs, and k.
118 verbosity : int (optional)
119 If greater than 0, prints solve time and number of iterations.
120 Each GP is created and solved with verbosity one less than this, so
121 if greater than 1, prints solver name and time for each GP.
122 x0 : dict (optional)
123 Initial location to approximate signomials about.
124 reltol : float
125 Iteration ends when this is greater than the distance between two
126 consecutive solve's objective values.
127 iteration_limit : int
128 Maximum GP iterations allowed.
129 mutategp: boolean
130 Prescribes whether to mutate the previously generated GP
131 or to create a new GP with every solve.
132 **solveargs :
133 Passed to solver function.
135 Returns
136 -------
137 result : dict
138 A dictionary containing the translated solver result.
139 """
140 self.gps, self.solver_outs, self._results = [], [], []
141 starttime = time()
142 if verbosity > 0:
143 print("Starting a sequence of GP solves")
144 print(f" for {len(self.sgpvks)} free variables")
145 print(f" in {len(self.sgpconstraints)} locally-GP constraints")
146 print(f" and for {len(self._gp.varlocs)} free variables")
147 print(f" in {len(self._gp.k)} posynomial inequalities.")
148 prevcost, cost, rel_improvement = None, None, None
149 while rel_improvement is None or rel_improvement > reltol:
150 prevcost = cost
151 if len(self.gps) > iteration_limit:
152 raise Infeasible(
153 f"Unsolved after {len(self.gps)} iterations. Check "
154 "`m.program.results`; if they're converging, try "
155 "`.localsolve(..., iteration_limit=NEWLIMIT)`.")
156 gp = self.gp(x0, cleanx0=len(self.gps) >= 1) # clean the first x0
157 self.gps.append(gp) # NOTE: SIDE EFFECTS
158 if verbosity > 1:
159 print(f"\nGP Solve {len(self.gps)}")
160 if verbosity > 2:
161 print("===============")
162 solver_out = gp.solve(solver, verbosity=verbosity-1,
163 gen_result=False, **solveargs)
164 self.solver_outs.append(solver_out)
165 cost = float(solver_out["objective"])
166 x0 = dict(zip(gp.varlocs, np.exp(solver_out["primal"])))
167 if verbosity > 2:
168 result = gp.generate_result(solver_out, verbosity=verbosity-3)
169 self._results.append(result)
170 vartable = result.table(self.sgpvks, tables=["freevariables"])
171 vartable = "\n" + vartable.replace("Free", "SGP", 1)
172 print(vartable)
173 elif verbosity > 1:
174 print(f"Solved cost was {cost:.4g}.")
175 if prevcost is None:
176 continue
177 rel_improvement = (prevcost - cost)/(prevcost + cost)
178 if cost/prevcost >= 1 + 10*EPS:
179 pywarnings.warn(
180 "SGP not convergent: Cost rose by "
181 f"{100*(cost - prevcost)/prevcost:.2g}% "
182 f"({prevcost:.6g} to {cost:.6g}) on GP solve "
183 f"{len(self.gps)}. Details can be found in "
184 "`m.program.results` or by solving at a higher "
185 "verbosity. Note convergence is not guaranteed for "
186 "models with SignomialEqualities.")
187 rel_improvement = cost = None
188 # solved successfully!
189 self.result = gp.generate_result(solver_out, verbosity=verbosity-3)
190 self.result["soltime"] = time() - starttime
191 if verbosity > 1:
192 print()
193 if verbosity > 0:
194 print(f"Solving took {self.result['soltime']:.3g} seconds and "
195 f"{len(self.gps)} GP solves.")
196 if hasattr(self.slack, "key"):
197 initsolwarning(self.result, "Slack Non-GP Constraints")
198 excess_slack = self.result["variables"][self.slack.key] - 1 # pylint: disable=no-member
199 if excess_slack > EPS:
200 msg = ("Final PCCP solution let non-GP constraints slacken by"
201 f" {100*excess_slack:.2g}%.")
202 expl = (msg +
203 " Calling .localsolve(pccp_penalty=...) with a higher"
204 f" `pccp_penalty` (it was {self.pccp_penalty:.3g} this"
205 " time) will reduce slack if the model is solvable"
206 "with less. To verify that the slack is needed, "
207 "generate an SGP with `use_pccp=False` and start it "
208 "from this model's solution: e.g. `m.localsolve("
209 "use_pccp=False, x0=m.solution[\"variables\"])`.")
210 if err_on_relax:
211 raise Infeasible(expl)
212 appendsolwarning(msg, (1 + excess_slack), self.result,
213 "Slack Non-GP Constraints")
214 if verbosity > -1:
215 print(expl)
216 self.result["cost function"] = self.cost
217 del self.result["freevariables"][self.slack.key] # pylint: disable=no-member
218 del self.result["variables"][self.slack.key] # pylint: disable=no-member
219 if "sensitivities" in self.result: # not true for MIGP
220 del self.result["sensitivities"]["variables"][self.slack.key] # pylint: disable=no-member
221 del self.result["sensitivities"]["variablerisk"][self.slack.key] # pylint: disable=no-member
222 slcon = self.gpconstraints[0]
223 slconsenss = self.result["sensitivities"]["constraints"][slcon]
224 del self.result["sensitivities"]["constraints"][slcon]
225 # TODO: create constraint in RelaxPCCP namespace
226 self.result["sensitivities"]["models"][""] -= slconsenss
227 if not self.result["sensitivities"]["models"][""]:
228 del self.result["sensitivities"]["models"][""]
229 return self.result
231 @property
232 def results(self):
233 "Creates and caches results from the raw solver_outs"
234 if not self._results:
235 self._results = [gp.generate_result(s_o, dual_check=False)
236 for gp, s_o in zip(self.gps, self.solver_outs)]
237 return self._results
239 def gp(self, x0=None, *, cleanx0=False):
240 "Update self._gp for x0 and return it."
241 if not x0:
242 return self._gp # return last generated
243 if not cleanx0:
244 cleanedx0 = KeyDict()
245 cleanedx0.vks = self._gp.x0.vks
246 cleanedx0.update(x0)
247 x0 = cleanedx0
248 self._gp.x0.update({vk: x0[vk] for vk in self.sgpvks if vk in x0})
249 p_idx = 0
250 for sgpc in self.sgpconstraints:
251 for hmaplt1 in sgpc.as_gpconstr(self._gp.x0).as_hmapslt1({}):
252 approxc = self.approxconstraints[p_idx]
253 approxc.left = self.slack
254 approxc.right.hmap = hmaplt1
255 approxc.unsubbed = [Posynomial(hmaplt1)/self.slack]
256 p_idx += 1 # p_idx=0 is the cost; sp constraints are after it
257 hmap, = approxc.as_hmapslt1(self._gp.substitutions)
258 self._gp.hmaps[p_idx] = hmap
259 m_idx = self._gp.m_idxs[p_idx].start
260 a_idxs = list(self.a_idxs[p_idx]) # A's entries we can modify
261 for i, (exp, c) in enumerate(hmap.items()):
262 self._gp.exps[m_idx + i] = exp
263 self._gp.cs[m_idx + i] = c
264 for var, x in exp.items():
265 try: # modify a particular A entry
266 row_idx = a_idxs.pop()
267 self._gp.A.row[row_idx] = m_idx + i
268 self._gp.A.col[row_idx] = self._gp.varidxs[var]
269 self._gp.A.data[row_idx] = x
270 except IndexError: # numbers of exps increased
271 self.a_idxs[p_idx].append(len(self._gp.A.row))
272 self._gp.A.row.append(m_idx + i)
273 self._gp.A.col.append(self._gp.varidxs[var])
274 self._gp.A.data.append(x)
275 for row_idx in a_idxs: # number of exps decreased
276 self._gp.A.row[row_idx] = 0 # zero out this entry
277 self._gp.A.col[row_idx] = 0
278 self._gp.A.data[row_idx] = 0
279 return self._gp