Coverage for gpkit/solvers/mosek_conif.py: 100%
142 statements
« prev ^ index » next coverage.py v7.4.0, created at 2024-01-05 22:33 -0500
« prev ^ index » next coverage.py v7.4.0, created at 2024-01-05 22:33 -0500
1"""Implements the GPkit interface to MOSEK (version >= 9)
2 through python-based Optimizer API"""
3import numpy as np
4import mosek # pylint: disable=import-error
5from ..exceptions import (UnknownInfeasible, InvalidLicense,
6 PrimalInfeasible, DualInfeasible)
8def optimize(*, c, A, k, p_idxs, **kwargs):
9 # pylint: disable=too-many-locals,too-many-statements,too-many-branches
10 # pylint: disable=invalid-name
11 """
12 Definitions
13 -----------
14 "[a,b] array of floats" indicates array-like data with shape [a,b]
15 n is the number of monomials in the gp
16 m is the number of variables in the gp
17 p is the number of posynomial constraints in the gp
19 Arguments
20 ---------
21 c : floats array of shape n
22 Coefficients of each monomial
23 A : gpkit.small_classes.CootMatrix, of shape (n, m)
24 Exponents of the various free variables for each monomial.
25 k : ints array of shape p+1
26 k[0] is the number of monomials (rows of A) present in the objective
27 k[1:] is the number of monomials (rows of A) present in each constraint
28 p_idxs : ints array of shape n.
29 sel = p_idxs == i selects rows of A
30 and entries of c of the i-th posynomial
31 fi(x) = c[sel] @ exp(A[sel,:] @ x).
32 The 0-th posynomial gives the objective function, and the
33 remaining posynomials should be constrained to be <= 1.
35 Returns
36 -------
37 dict
38 Contains the following keys
39 "status": string
40 "optimal", "infeasible", "unbounded", or "unknown".
41 "primal" np.ndarray or None
42 The values of the ``m`` primal variables.
43 "la": np.ndarray or None
44 The dual variables to the ``p`` posynomial constraints, when
45 those constraints are represented in log-sum-exp ("LSE") form.
46 """
47 #
48 # Initial transformations of problem data.
49 #
50 # separate monomial constraints (call them "lin"),
51 # from those which require
52 # an LSE representation (call those "lse").
53 #
54 # NOTE: the epigraph of the objective function always gets an "lse"
55 # representation, even if the objective is a monomial.
56 #
57 log_c = np.log(np.array(c))
58 lse_posys = [0]
59 lin_posys = []
60 for i, val in enumerate(k[1:]):
61 if val > 1:
62 lse_posys.append(i+1)
63 else:
64 lin_posys.append(i+1)
65 if lin_posys:
66 A = A.tocsr()
67 lin_posys_set = frozenset(lin_posys)
68 lin_idxs = [i for i, p in enumerate(p_idxs) if p in lin_posys_set]
69 lse_idxs = np.ones(A.shape[0], dtype=bool)
70 lse_idxs[lin_idxs] = False
71 A_lse = A[lse_idxs, :].tocoo()
72 log_c_lse = log_c[lse_idxs]
73 A_lin = A[lin_idxs, :].tocoo()
74 log_c_lin = log_c[lin_idxs]
75 else:
76 log_c_lin = None # A_lin won't be referenced later,
77 A_lse = A # so no need to define it.
78 log_c_lse = log_c
79 k_lse = [k[i] for i in lse_posys]
80 n_lse = sum(k_lse)
81 p_lse = len(k_lse)
82 lse_p_idx = []
83 for i, ki in enumerate(k_lse):
84 lse_p_idx.extend([i] * ki)
85 lse_p_idx = np.array(lse_p_idx)
86 #
87 # Create MOSEK task. Add variables, bound constraints, and conic
88 # constraints.
89 #
90 # Say that MOSEK's optimization variable is a block vector,
91 # [x;t;z], where ...
92 # x is the user-defined primal variable (length m),
93 # t is an aux variable for exponential cones (length 3 * n_lse),
94 # z is an epigraph variable for LSE terms (length p_lse).
95 #
96 # The variable z[0] is special,
97 # because it's the epigraph of the objective function
98 # in LSE form. The sign of this variable is not constrained.
99 #
100 # The variables z[1:] are epigraph terms for "log",
101 # in constraints that naturally write as
102 # LSE(Ai @ x + log_ci) <= 0. These variables need to be <= 0.
103 #
104 # The main constraints on (x, t, z) are described
105 # in next comment block.
106 #
107 env = mosek.Env()
108 task = env.Task(0, 0)
109 m = A.shape[1]
110 msk_nvars = m + 3 * n_lse + p_lse
111 task.appendvars(msk_nvars)
112 # "x" is free
113 task.putvarboundlist(np.arange(m), [mosek.boundkey.fr] * m,
114 np.zeros(m), np.zeros(m))
115 # t[3 * i + i] == 1, other components are free.
116 bound_types = [mosek.boundkey.fr, mosek.boundkey.fx, mosek.boundkey.fr]
117 task.putvarboundlist(np.arange(m, m + 3*n_lse), bound_types * n_lse,
118 np.ones(3*n_lse), np.ones(3*n_lse))
119 # z[0] is free; z[1:] <= 0.
120 bound_types = [mosek.boundkey.fr] + [mosek.boundkey.up] * (p_lse - 1)
121 task.putvarboundlist(np.arange(m + 3*n_lse, msk_nvars), bound_types,
122 np.zeros(p_lse), np.zeros(p_lse))
123 # t[3*i], t[3*i + 1], t[3*i + 2] belongs to the exponential cone
124 task.appendconesseq([mosek.conetype.pexp] * n_lse, [0.0] * n_lse,
125 [3] * n_lse, m)
126 #
127 # Exponential cone affine constraints (other than t[3*i + 1] == 1).
128 #
129 # For each i in {0, ..., n_lse - 1}, we need
130 # t[3*i + 2] == A_lse[i, :] @ x + log_c_lse[i] - z[lse_p_idx[i]].
131 #
132 # For each j from {0, ..., p_lse - 1}, the "t" should also satisfy
133 # sum(t[3*i] for i where i == lse_p_idx[j]) <= 1.
134 #
135 # When combined with bound constraints ``t[3*i + 1] == 1``, the
136 # above constraints imply
137 # LSE(A_lse[sel, :] @ x + log_c_lse[sel]) <= z[i]
138 # for ``sel = lse_p_idx == i``.
139 #
140 task.appendcons(n_lse + p_lse)
141 # Linear equations between (x,t,z).
142 # start with coefficients on "x"
143 rows = list(A_lse.row)
144 cols = list(A_lse.col)
145 vals = list(A_lse.data)
146 # add coefficients on "t"
147 rows += list(range(n_lse))
148 cols += (m + 3*np.arange(n_lse) + 2).tolist()
149 vals += [-1.0] * n_lse
150 # add coefficients on "z"
151 rows += list(range(n_lse))
152 cols += [m + 3*n_lse + lse_p_idx[i] for i in range(n_lse)]
153 vals += [-1.0] * n_lse
154 task.putaijlist(rows, cols, vals)
155 cur_con_idx = n_lse
156 # Linear inequalities on certain sums of "t".
157 rows, cols, vals = [], [], []
158 for i in range(p_lse):
159 sels = np.nonzero(lse_p_idx == i)[0]
160 rows.extend([cur_con_idx] * sels.size)
161 cols.extend(m + 3 * sels)
162 vals.extend([1] * sels.size)
163 cur_con_idx += 1
164 task.putaijlist(rows, cols, vals)
165 # Build the right-hand-sides of the [in]equality constraints
166 type_constraint = [mosek.boundkey.fx] * n_lse + [mosek.boundkey.up] * p_lse
167 h = np.concatenate([-log_c_lse, np.ones(p_lse)])
168 task.putconboundlist(np.arange(h.size, dtype=int), type_constraint, h, h)
169 #
170 # Affine constraints, not needing the exponential cone
171 #
172 # Require A_lin @ x <= -log_c_lin.
173 #
174 if log_c_lin is not None:
175 task.appendcons(log_c_lin.size)
176 rows = cur_con_idx + np.array(A_lin.row)
177 task.putaijlist(rows, A_lin.col, A_lin.data)
178 type_constraint = [mosek.boundkey.up] * log_c_lin.size
179 con_indices = np.arange(cur_con_idx, cur_con_idx + log_c_lin.size)
180 h = -log_c_lin #pylint: disable=invalid-unary-operand-type
181 task.putconboundlist(con_indices, type_constraint, h, h)
182 cur_con_idx += log_c_lin.size
183 #
184 # Constrain choice variables to discrete choices
185 #
186 choicevaridxs = kwargs.get("choicevaridxs", {})
187 if choicevaridxs:
188 n_choicevars = 0
189 for var, idx in choicevaridxs.items():
190 choices = sorted(var.choices)
191 m_choices = len(choices) - 1 # the first option is the default
192 choiceidxs = np.arange(msk_nvars + n_choicevars,
193 msk_nvars + n_choicevars + m_choices)
194 n_choicevars += m_choices
195 task.appendvars(m_choices)
196 task.putvartypelist(choiceidxs,
197 [mosek.variabletype.type_int]*m_choices)
198 task.putvarboundlist(choiceidxs, [mosek.boundkey.ra]*m_choices,
199 np.zeros(m_choices), np.ones(m_choices))
200 task.appendcons(m_choices)
201 for i in range(m_choices - 1):
202 # each larger choice requires those before it
203 task.putarow(cur_con_idx + i, choiceidxs[i:i+2], [1.0, -1.0])
204 task.putconbound(cur_con_idx + i, mosek.boundkey.lo, 0.0, 0.0)
205 cur_con_idx += 1
206 base = np.log(choices[0])
207 logdiffs = np.diff(np.log(choices)).tolist()
208 task.putarow(cur_con_idx, choiceidxs.tolist() + [idx],
209 logdiffs + [-1]) # choices are summed
210 task.putconbound(cur_con_idx, mosek.boundkey.fx, -base, -base)
211 cur_con_idx += 1
212 #
213 # Set the objective function
214 #
215 task.putclist([int(m + 3*n_lse)], [1])
216 task.putobjsense(mosek.objsense.minimize)
217 #
218 # Set solver parameters, and call .solve().
219 #
220 verbose = kwargs.get("verbose", True)
221 if verbose:
222 def streamprinter(text):
223 "Stream printer for output from mosek."
224 print(text)
226 env.set_Stream(mosek.streamtype.log, streamprinter)
227 task.set_Stream(mosek.streamtype.log, streamprinter)
228 task.putintparam(mosek.iparam.infeas_report_auto, mosek.onoffkey.on)
229 task.putintparam(mosek.iparam.log_presolve, 0)
231 try:
232 task.optimize()
233 except mosek.Error as e: # pragma: no cover
234 if e.errno in [mosek.rescode.err_missing_license_file,
235 mosek.rescode.err_license_version,
236 mosek.rescode.err_license_expired]:
237 raise InvalidLicense() from e
238 raise e
240 if verbose:
241 task.solutionsummary(mosek.streamtype.msg)
242 #
243 # Recover the solution
244 #
245 if choicevaridxs:
246 sol = mosek.soltype.itg
247 optimal = mosek.solsta.integer_optimal
248 else:
249 sol = mosek.soltype.itr
250 optimal = mosek.solsta.optimal
251 msk_solsta = task.getsolsta(sol)
252 if msk_solsta == mosek.solsta.prim_infeas_cer:
253 raise PrimalInfeasible()
254 if msk_solsta == mosek.solsta.dual_infeas_cer:
255 raise DualInfeasible()
256 if msk_solsta != optimal: # pragma: no cover
257 raise UnknownInfeasible("solution status: ", msk_solsta)
259 # recover primal variables
260 x = [0.] * m
261 task.getxxslice(sol, 0, m, x)
262 # recover binary variables
263 # xbin = [0.] * (n_choicevars)
264 # task.getxxslice(sol, msk_nvars, msk_nvars + n_choicevars, xbin)
265 # wrap things up in a dictionary
266 solution = {"status": "optimal", "primal": np.array(x),
267 "objective": np.exp(task.getprimalobj(sol))}
268 # recover dual variables for log-sum-exp epigraph constraints
269 # (skip epigraph of the objective function).
270 if choicevaridxs: # no dual solution
271 solution["la"] = []
272 solution["nu"] = []
273 else:
274 z_duals = [0.] * (p_lse - 1)
275 task.getsuxslice(mosek.soltype.itr, m + 3*n_lse + 1, msk_nvars, z_duals)
276 z_duals = np.array(z_duals)
277 z_duals[z_duals < 0] = 0
278 # recover dual variables for the remaining user-provided constraints
279 if log_c_lin is None:
280 solution["la"] = z_duals
281 else:
282 aff_duals = [0.] * log_c_lin.size
283 task.getsucslice(mosek.soltype.itr, n_lse + p_lse, cur_con_idx,
284 aff_duals)
285 aff_duals = np.array(aff_duals)
286 aff_duals[aff_duals < 0] = 0
287 # merge z_duals with aff_duals
288 merged_duals = np.zeros(len(k))
289 merged_duals[lse_posys[1:]] = z_duals # skipping the cost
290 merged_duals[lin_posys] = aff_duals
291 solution["la"] = merged_duals[1:]
293 task.__exit__(None, None, None)
294 env.__exit__(None, None, None)
295 return solution