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