Coverage for gpkit/solvers/mosek_conif.py: 100%

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

142 statements  

1"""Implements the GPkit interface to MOSEK (version >= 9) 

2 through python-based Optimizer API""" 

3import numpy as np 

4import mosek 

5from ..exceptions import (UnknownInfeasible, InvalidLicense, 

6 PrimalInfeasible, DualInfeasible) 

7 

8def optimize(*, c, A, k, p_idxs, **kwargs): 

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

10 """ 

11 Definitions 

12 ----------- 

13 "[a,b] array of floats" indicates array-like data with shape [a,b] 

14 n is the number of monomials in the gp 

15 m is the number of variables in the gp 

16 p is the number of posynomial constraints in the gp 

17 

18 Arguments 

19 --------- 

20 c : floats array of shape n 

21 Coefficients of each monomial 

22 A : gpkit.small_classes.CootMatrix, of shape (n, m) 

23 Exponents of the various free variables for each monomial. 

24 k : ints array of shape p+1 

25 k[0] is the number of monomials (rows of A) present in the objective 

26 k[1:] is the number of monomials (rows of A) present in each constraint 

27 p_idxs : ints array of shape n. 

28 sel = p_idxs == i selects rows of A 

29 and entries of c of the i-th posynomial 

30 fi(x) = c[sel] @ exp(A[sel,:] @ x). 

31 The 0-th posynomial gives the objective function, and the 

32 remaining posynomials should be constrained to be <= 1. 

33 

34 Returns 

35 ------- 

36 dict 

37 Contains the following keys 

38 "status": string 

39 "optimal", "infeasible", "unbounded", or "unknown". 

40 "primal" np.ndarray or None 

41 The values of the ``m`` primal variables. 

42 "la": np.ndarray or None 

43 The dual variables to the ``p`` posynomial constraints, when 

44 those constraints are represented in log-sum-exp ("LSE") form. 

45 """ 

46 # 

47 # Initial transformations of problem data. 

48 # 

49 # separate monomial constraints (call them "lin"), 

50 # from those which require 

51 # an LSE representation (call those "lse"). 

52 # 

53 # NOTE: the epigraph of the objective function always gets an "lse" 

54 # representation, even if the objective is a monomial. 

55 # 

56 log_c = np.log(np.array(c)) 

57 lse_posys = [0] 

58 lin_posys = [] 

59 for i, val in enumerate(k[1:]): 

60 if val > 1: 

61 lse_posys.append(i+1) 

62 else: 

63 lin_posys.append(i+1) 

64 if lin_posys: 

65 A = A.tocsr() 

66 lin_posys_set = frozenset(lin_posys) 

67 lin_idxs = [i for i, p in enumerate(p_idxs) if p in lin_posys_set] 

68 lse_idxs = np.ones(A.shape[0], dtype=bool) 

69 lse_idxs[lin_idxs] = False 

70 A_lse = A[lse_idxs, :].tocoo() 

71 log_c_lse = log_c[lse_idxs] 

72 A_lin = A[lin_idxs, :].tocoo() 

73 log_c_lin = log_c[lin_idxs] 

74 else: 

75 log_c_lin = None # A_lin won't be referenced later, 

76 A_lse = A # so no need to define it. 

77 log_c_lse = log_c 

78 k_lse = [k[i] for i in lse_posys] 

79 n_lse = sum(k_lse) 

80 p_lse = len(k_lse) 

81 lse_p_idx = [] 

82 for i, ki in enumerate(k_lse): 

83 lse_p_idx.extend([i] * ki) 

84 lse_p_idx = np.array(lse_p_idx) 

85 # 

86 # Create MOSEK task. Add variables, bound constraints, and conic 

87 # constraints. 

88 # 

89 # Say that MOSEK's optimization variable is a block vector, 

90 # [x;t;z], where ... 

91 # x is the user-defined primal variable (length m), 

92 # t is an aux variable for exponential cones (length 3 * n_lse), 

93 # z is an epigraph variable for LSE terms (length p_lse). 

94 # 

95 # The variable z[0] is special, 

96 # because it's the epigraph of the objective function 

97 # in LSE form. The sign of this variable is not constrained. 

98 # 

99 # The variables z[1:] are epigraph terms for "log", 

100 # in constraints that naturally write as 

101 # LSE(Ai @ x + log_ci) <= 0. These variables need to be <= 0. 

102 # 

103 # The main constraints on (x, t, z) are described 

104 # in next comment block. 

105 # 

106 env = mosek.Env() 

107 task = env.Task(0, 0) 

108 m = A.shape[1] 

109 msk_nvars = m + 3 * n_lse + p_lse 

110 task.appendvars(msk_nvars) 

111 # "x" is free 

112 task.putvarboundlist(np.arange(m), [mosek.boundkey.fr] * m, 

113 np.zeros(m), np.zeros(m)) 

114 # t[3 * i + i] == 1, other components are free. 

115 bound_types = [mosek.boundkey.fr, mosek.boundkey.fx, mosek.boundkey.fr] 

116 task.putvarboundlist(np.arange(m, m + 3*n_lse), bound_types * n_lse, 

117 np.ones(3*n_lse), np.ones(3*n_lse)) 

118 # z[0] is free; z[1:] <= 0. 

119 bound_types = [mosek.boundkey.fr] + [mosek.boundkey.up] * (p_lse - 1) 

120 task.putvarboundlist(np.arange(m + 3*n_lse, msk_nvars), bound_types, 

121 np.zeros(p_lse), np.zeros(p_lse)) 

122 # t[3*i], t[3*i + 1], t[3*i + 2] belongs to the exponential cone 

123 task.appendconesseq([mosek.conetype.pexp] * n_lse, [0.0] * n_lse, 

124 [3] * n_lse, m) 

125 # 

126 # Exponential cone affine constraints (other than t[3*i + 1] == 1). 

127 # 

128 # For each i in {0, ..., n_lse - 1}, we need 

129 # t[3*i + 2] == A_lse[i, :] @ x + log_c_lse[i] - z[lse_p_idx[i]]. 

130 # 

131 # For each j from {0, ..., p_lse - 1}, the "t" should also satisfy 

132 # sum(t[3*i] for i where i == lse_p_idx[j]) <= 1. 

133 # 

134 # When combined with bound constraints ``t[3*i + 1] == 1``, the 

135 # above constraints imply 

136 # LSE(A_lse[sel, :] @ x + log_c_lse[sel]) <= z[i] 

137 # for ``sel = lse_p_idx == i``. 

138 # 

139 task.appendcons(n_lse + p_lse) 

140 # Linear equations between (x,t,z). 

141 # start with coefficients on "x" 

142 rows = list(A_lse.row) 

143 cols = list(A_lse.col) 

144 vals = list(A_lse.data) 

145 # add coefficients on "t" 

146 rows += list(range(n_lse)) 

147 cols += (m + 3*np.arange(n_lse) + 2).tolist() 

148 vals += [-1.0] * n_lse 

149 # add coefficients on "z" 

150 rows += list(range(n_lse)) 

151 cols += [m + 3*n_lse + lse_p_idx[i] for i in range(n_lse)] 

152 vals += [-1.0] * n_lse 

153 task.putaijlist(rows, cols, vals) 

154 cur_con_idx = n_lse 

155 # Linear inequalities on certain sums of "t". 

156 rows, cols, vals = [], [], [] 

157 for i in range(p_lse): 

158 sels = np.nonzero(lse_p_idx == i)[0] 

159 rows.extend([cur_con_idx] * sels.size) 

160 cols.extend(m + 3 * sels) 

161 vals.extend([1] * sels.size) 

162 cur_con_idx += 1 

163 task.putaijlist(rows, cols, vals) 

164 # Build the right-hand-sides of the [in]equality constraints 

165 type_constraint = [mosek.boundkey.fx] * n_lse + [mosek.boundkey.up] * p_lse 

166 h = np.concatenate([-log_c_lse, np.ones(p_lse)]) 

167 task.putconboundlist(np.arange(h.size, dtype=int), type_constraint, h, h) 

168 # 

169 # Affine constraints, not needing the exponential cone 

170 # 

171 # Require A_lin @ x <= -log_c_lin. 

172 # 

173 if log_c_lin is not None: 

174 task.appendcons(log_c_lin.size) 

175 rows = cur_con_idx + np.array(A_lin.row) 

176 task.putaijlist(rows, A_lin.col, A_lin.data) 

177 type_constraint = [mosek.boundkey.up] * log_c_lin.size 

178 con_indices = np.arange(cur_con_idx, cur_con_idx + log_c_lin.size) 

179 h = -log_c_lin #pylint: disable=invalid-unary-operand-type 

180 task.putconboundlist(con_indices, type_constraint, h, h) 

181 cur_con_idx += log_c_lin.size 

182 # 

183 # Constrain choice variables to discrete choices 

184 # 

185 choicevaridxs = kwargs.get("choicevaridxs", {}) 

186 if choicevaridxs: 

187 n_choicevars = 0 

188 for var, idx in choicevaridxs.items(): 

189 choices = sorted(var.choices) 

190 m_choices = len(choices) - 1 # the first option is the default 

191 choiceidxs = np.arange(msk_nvars + n_choicevars, 

192 msk_nvars + n_choicevars + m_choices) 

193 n_choicevars += m_choices 

194 task.appendvars(m_choices) 

195 task.putvartypelist(choiceidxs, 

196 [mosek.variabletype.type_int]*m_choices) 

197 task.putvarboundlist(choiceidxs, [mosek.boundkey.ra]*m_choices, 

198 np.zeros(m_choices), np.ones(m_choices)) 

199 task.appendcons(m_choices) 

200 for i in range(m_choices - 1): 

201 # each larger choice requires those before it 

202 task.putarow(cur_con_idx + i, choiceidxs[i:i+2], [1.0, -1.0]) 

203 task.putconbound(cur_con_idx + i, mosek.boundkey.lo, 0.0, 0.0) 

204 cur_con_idx += 1 

205 base = np.log(choices[0]) 

206 logdiffs = np.diff(np.log(choices)).tolist() 

207 task.putarow(cur_con_idx, choiceidxs.tolist() + [idx], 

208 logdiffs + [-1]) # choices are summed 

209 task.putconbound(cur_con_idx, mosek.boundkey.fx, -base, -base) 

210 cur_con_idx += 1 

211 # 

212 # Set the objective function 

213 # 

214 task.putclist([int(m + 3*n_lse)], [1]) 

215 task.putobjsense(mosek.objsense.minimize) 

216 # 

217 # Set solver parameters, and call .solve(). 

218 # 

219 verbose = kwargs.get("verbose", True) 

220 if verbose: 

221 def streamprinter(text): 

222 "Stream printer for output from mosek." 

223 print(text) 

224 

225 env.set_Stream(mosek.streamtype.log, streamprinter) 

226 task.set_Stream(mosek.streamtype.log, streamprinter) 

227 task.putintparam(mosek.iparam.infeas_report_auto, mosek.onoffkey.on) 

228 task.putintparam(mosek.iparam.log_presolve, 0) 

229 

230 try: 

231 task.optimize() 

232 except mosek.Error as e: # pragma: no cover 

233 if e.errno in [mosek.rescode.err_missing_license_file, 

234 mosek.rescode.err_license_version, 

235 mosek.rescode.err_license_expired]: 

236 raise InvalidLicense() from e 

237 raise e 

238 

239 if verbose: 

240 task.solutionsummary(mosek.streamtype.msg) 

241 # 

242 # Recover the solution 

243 # 

244 if choicevaridxs: 

245 sol = mosek.soltype.itg 

246 optimal = mosek.solsta.integer_optimal 

247 else: 

248 sol = mosek.soltype.itr 

249 optimal = mosek.solsta.optimal 

250 msk_solsta = task.getsolsta(sol) 

251 if msk_solsta == mosek.solsta.prim_infeas_cer: 

252 raise PrimalInfeasible() 

253 if msk_solsta == mosek.solsta.dual_infeas_cer: 

254 raise DualInfeasible() 

255 if msk_solsta != optimal: # pragma: no cover 

256 raise UnknownInfeasible("solution status: ", msk_solsta) 

257 

258 # recover primal variables 

259 x = [0.] * m 

260 task.getxxslice(sol, 0, m, x) 

261 # recover binary variables 

262 # xbin = [0.] * (n_choicevars) 

263 # task.getxxslice(sol, msk_nvars, msk_nvars + n_choicevars, xbin) 

264 # wrap things up in a dictionary 

265 solution = {"status": "optimal", "primal": np.array(x), 

266 "objective": np.exp(task.getprimalobj(sol))} 

267 # recover dual variables for log-sum-exp epigraph constraints 

268 # (skip epigraph of the objective function). 

269 if choicevaridxs: # no dual solution 

270 solution["la"] = [] 

271 solution["nu"] = [] 

272 else: 

273 z_duals = [0.] * (p_lse - 1) 

274 task.getsuxslice(mosek.soltype.itr, m + 3*n_lse + 1, msk_nvars, z_duals) 

275 z_duals = np.array(z_duals) 

276 z_duals[z_duals < 0] = 0 

277 # recover dual variables for the remaining user-provided constraints 

278 if log_c_lin is None: 

279 solution["la"] = z_duals 

280 else: 

281 aff_duals = [0.] * log_c_lin.size 

282 task.getsucslice(mosek.soltype.itr, n_lse + p_lse, cur_con_idx, 

283 aff_duals) 

284 aff_duals = np.array(aff_duals) 

285 aff_duals[aff_duals < 0] = 0 

286 # merge z_duals with aff_duals 

287 merged_duals = np.zeros(len(k)) 

288 merged_duals[lse_posys[1:]] = z_duals # skipping the cost 

289 merged_duals[lin_posys] = aff_duals 

290 solution["la"] = merged_duals[1:] 

291 

292 task.__exit__(None, None, None) 

293 env.__exit__(None, None, None) 

294 return solution