Hide keyboard shortcuts

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""" 

2This module implements the Sequential Least Squares Programming optimization 

3algorithm (SLSQP), originally developed by Dieter Kraft. 

4See http://www.netlib.org/toms/733 

5 

6Functions 

7--------- 

8.. autosummary:: 

9 :toctree: generated/ 

10 

11 approx_jacobian 

12 fmin_slsqp 

13 

14""" 

15 

16__all__ = ['approx_jacobian', 'fmin_slsqp'] 

17 

18import numpy as np 

19from scipy.optimize._slsqp import slsqp 

20from numpy import (zeros, array, linalg, append, asfarray, concatenate, finfo, 

21 sqrt, vstack, exp, inf, isfinite, atleast_1d) 

22from .optimize import (OptimizeResult, _check_unknown_options, 

23 _prepare_scalar_function) 

24from ._numdiff import approx_derivative 

25from ._constraints import old_bound_to_new 

26 

27 

28__docformat__ = "restructuredtext en" 

29 

30_epsilon = sqrt(finfo(float).eps) 

31 

32 

33def approx_jacobian(x, func, epsilon, *args): 

34 """ 

35 Approximate the Jacobian matrix of a callable function. 

36 

37 Parameters 

38 ---------- 

39 x : array_like 

40 The state vector at which to compute the Jacobian matrix. 

41 func : callable f(x,*args) 

42 The vector-valued function. 

43 epsilon : float 

44 The perturbation used to determine the partial derivatives. 

45 args : sequence 

46 Additional arguments passed to func. 

47 

48 Returns 

49 ------- 

50 An array of dimensions ``(lenf, lenx)`` where ``lenf`` is the length 

51 of the outputs of `func`, and ``lenx`` is the number of elements in 

52 `x`. 

53 

54 Notes 

55 ----- 

56 The approximation is done using forward differences. 

57 

58 """ 

59 # approx_derivative returns (m, n) == (lenf, lenx) 

60 jac = approx_derivative(func, x, method='2-point', abs_step=epsilon, 

61 args=args) 

62 # if func returns a scalar jac.shape will be (lenx,). Make sure 

63 # it's at least a 2D array. 

64 return np.atleast_2d(jac) 

65 

66 

67def fmin_slsqp(func, x0, eqcons=(), f_eqcons=None, ieqcons=(), f_ieqcons=None, 

68 bounds=(), fprime=None, fprime_eqcons=None, 

69 fprime_ieqcons=None, args=(), iter=100, acc=1.0E-6, 

70 iprint=1, disp=None, full_output=0, epsilon=_epsilon, 

71 callback=None): 

72 """ 

73 Minimize a function using Sequential Least Squares Programming 

74 

75 Python interface function for the SLSQP Optimization subroutine 

76 originally implemented by Dieter Kraft. 

77 

78 Parameters 

79 ---------- 

80 func : callable f(x,*args) 

81 Objective function. Must return a scalar. 

82 x0 : 1-D ndarray of float 

83 Initial guess for the independent variable(s). 

84 eqcons : list, optional 

85 A list of functions of length n such that 

86 eqcons[j](x,*args) == 0.0 in a successfully optimized 

87 problem. 

88 f_eqcons : callable f(x,*args), optional 

89 Returns a 1-D array in which each element must equal 0.0 in a 

90 successfully optimized problem. If f_eqcons is specified, 

91 eqcons is ignored. 

92 ieqcons : list, optional 

93 A list of functions of length n such that 

94 ieqcons[j](x,*args) >= 0.0 in a successfully optimized 

95 problem. 

96 f_ieqcons : callable f(x,*args), optional 

97 Returns a 1-D ndarray in which each element must be greater or 

98 equal to 0.0 in a successfully optimized problem. If 

99 f_ieqcons is specified, ieqcons is ignored. 

100 bounds : list, optional 

101 A list of tuples specifying the lower and upper bound 

102 for each independent variable [(xl0, xu0),(xl1, xu1),...] 

103 Infinite values will be interpreted as large floating values. 

104 fprime : callable `f(x,*args)`, optional 

105 A function that evaluates the partial derivatives of func. 

106 fprime_eqcons : callable `f(x,*args)`, optional 

107 A function of the form `f(x, *args)` that returns the m by n 

108 array of equality constraint normals. If not provided, 

109 the normals will be approximated. The array returned by 

110 fprime_eqcons should be sized as ( len(eqcons), len(x0) ). 

111 fprime_ieqcons : callable `f(x,*args)`, optional 

112 A function of the form `f(x, *args)` that returns the m by n 

113 array of inequality constraint normals. If not provided, 

114 the normals will be approximated. The array returned by 

115 fprime_ieqcons should be sized as ( len(ieqcons), len(x0) ). 

116 args : sequence, optional 

117 Additional arguments passed to func and fprime. 

118 iter : int, optional 

119 The maximum number of iterations. 

120 acc : float, optional 

121 Requested accuracy. 

122 iprint : int, optional 

123 The verbosity of fmin_slsqp : 

124 

125 * iprint <= 0 : Silent operation 

126 * iprint == 1 : Print summary upon completion (default) 

127 * iprint >= 2 : Print status of each iterate and summary 

128 disp : int, optional 

129 Overrides the iprint interface (preferred). 

130 full_output : bool, optional 

131 If False, return only the minimizer of func (default). 

132 Otherwise, output final objective function and summary 

133 information. 

134 epsilon : float, optional 

135 The step size for finite-difference derivative estimates. 

136 callback : callable, optional 

137 Called after each iteration, as ``callback(x)``, where ``x`` is the 

138 current parameter vector. 

139 

140 Returns 

141 ------- 

142 out : ndarray of float 

143 The final minimizer of func. 

144 fx : ndarray of float, if full_output is true 

145 The final value of the objective function. 

146 its : int, if full_output is true 

147 The number of iterations. 

148 imode : int, if full_output is true 

149 The exit mode from the optimizer (see below). 

150 smode : string, if full_output is true 

151 Message describing the exit mode from the optimizer. 

152 

153 See also 

154 -------- 

155 minimize: Interface to minimization algorithms for multivariate 

156 functions. See the 'SLSQP' `method` in particular. 

157 

158 Notes 

159 ----- 

160 Exit modes are defined as follows :: 

161 

162 -1 : Gradient evaluation required (g & a) 

163 0 : Optimization terminated successfully 

164 1 : Function evaluation required (f & c) 

165 2 : More equality constraints than independent variables 

166 3 : More than 3*n iterations in LSQ subproblem 

167 4 : Inequality constraints incompatible 

168 5 : Singular matrix E in LSQ subproblem 

169 6 : Singular matrix C in LSQ subproblem 

170 7 : Rank-deficient equality constraint subproblem HFTI 

171 8 : Positive directional derivative for linesearch 

172 9 : Iteration limit reached 

173 

174 Examples 

175 -------- 

176 Examples are given :ref:`in the tutorial <tutorial-sqlsp>`. 

177 

178 """ 

179 if disp is not None: 

180 iprint = disp 

181 

182 opts = {'maxiter': iter, 

183 'ftol': acc, 

184 'iprint': iprint, 

185 'disp': iprint != 0, 

186 'eps': epsilon, 

187 'callback': callback} 

188 

189 # Build the constraints as a tuple of dictionaries 

190 cons = () 

191 # 1. constraints of the 1st kind (eqcons, ieqcons); no Jacobian; take 

192 # the same extra arguments as the objective function. 

193 cons += tuple({'type': 'eq', 'fun': c, 'args': args} for c in eqcons) 

194 cons += tuple({'type': 'ineq', 'fun': c, 'args': args} for c in ieqcons) 

195 # 2. constraints of the 2nd kind (f_eqcons, f_ieqcons) and their Jacobian 

196 # (fprime_eqcons, fprime_ieqcons); also take the same extra arguments 

197 # as the objective function. 

198 if f_eqcons: 

199 cons += ({'type': 'eq', 'fun': f_eqcons, 'jac': fprime_eqcons, 

200 'args': args}, ) 

201 if f_ieqcons: 

202 cons += ({'type': 'ineq', 'fun': f_ieqcons, 'jac': fprime_ieqcons, 

203 'args': args}, ) 

204 

205 res = _minimize_slsqp(func, x0, args, jac=fprime, bounds=bounds, 

206 constraints=cons, **opts) 

207 if full_output: 

208 return res['x'], res['fun'], res['nit'], res['status'], res['message'] 

209 else: 

210 return res['x'] 

211 

212 

213def _minimize_slsqp(func, x0, args=(), jac=None, bounds=None, 

214 constraints=(), 

215 maxiter=100, ftol=1.0E-6, iprint=1, disp=False, 

216 eps=_epsilon, callback=None, finite_diff_rel_step=None, 

217 **unknown_options): 

218 """ 

219 Minimize a scalar function of one or more variables using Sequential 

220 Least Squares Programming (SLSQP). 

221 

222 Options 

223 ------- 

224 ftol : float 

225 Precision goal for the value of f in the stopping criterion. 

226 eps : float 

227 Step size used for numerical approximation of the Jacobian. 

228 disp : bool 

229 Set to True to print convergence messages. If False, 

230 `verbosity` is ignored and set to 0. 

231 maxiter : int 

232 Maximum number of iterations. 

233 finite_diff_rel_step : None or array_like, optional 

234 If `jac in ['2-point', '3-point', 'cs']` the relative step size to 

235 use for numerical approximation of `jac`. The absolute step 

236 size is computed as ``h = rel_step * sign(x0) * max(1, abs(x0))``, 

237 possibly adjusted to fit into the bounds. For ``method='3-point'`` 

238 the sign of `h` is ignored. If None (default) then step is selected 

239 automatically. 

240 """ 

241 _check_unknown_options(unknown_options) 

242 iter = maxiter - 1 

243 acc = ftol 

244 epsilon = eps 

245 

246 if not disp: 

247 iprint = 0 

248 

249 # Constraints are triaged per type into a dictionary of tuples 

250 if isinstance(constraints, dict): 

251 constraints = (constraints, ) 

252 

253 cons = {'eq': (), 'ineq': ()} 

254 for ic, con in enumerate(constraints): 

255 # check type 

256 try: 

257 ctype = con['type'].lower() 

258 except KeyError: 

259 raise KeyError('Constraint %d has no type defined.' % ic) 

260 except TypeError: 

261 raise TypeError('Constraints must be defined using a ' 

262 'dictionary.') 

263 except AttributeError: 

264 raise TypeError("Constraint's type must be a string.") 

265 else: 

266 if ctype not in ['eq', 'ineq']: 

267 raise ValueError("Unknown constraint type '%s'." % con['type']) 

268 

269 # check function 

270 if 'fun' not in con: 

271 raise ValueError('Constraint %d has no function defined.' % ic) 

272 

273 # check Jacobian 

274 cjac = con.get('jac') 

275 if cjac is None: 

276 # approximate Jacobian function. The factory function is needed 

277 # to keep a reference to `fun`, see gh-4240. 

278 def cjac_factory(fun): 

279 def cjac(x, *args): 

280 if jac in ['2-point', '3-point', 'cs']: 

281 return approx_derivative(fun, x, method=jac, args=args, 

282 rel_step=finite_diff_rel_step) 

283 else: 

284 return approx_derivative(fun, x, method='2-point', 

285 abs_step=epsilon, args=args) 

286 

287 return cjac 

288 cjac = cjac_factory(con['fun']) 

289 

290 # update constraints' dictionary 

291 cons[ctype] += ({'fun': con['fun'], 

292 'jac': cjac, 

293 'args': con.get('args', ())}, ) 

294 

295 exit_modes = {-1: "Gradient evaluation required (g & a)", 

296 0: "Optimization terminated successfully", 

297 1: "Function evaluation required (f & c)", 

298 2: "More equality constraints than independent variables", 

299 3: "More than 3*n iterations in LSQ subproblem", 

300 4: "Inequality constraints incompatible", 

301 5: "Singular matrix E in LSQ subproblem", 

302 6: "Singular matrix C in LSQ subproblem", 

303 7: "Rank-deficient equality constraint subproblem HFTI", 

304 8: "Positive directional derivative for linesearch", 

305 9: "Iteration limit reached"} 

306 

307 # Transform x0 into an array. 

308 x = asfarray(x0).flatten() 

309 

310 # SLSQP is sent 'old-style' bounds, 'new-style' bounds are required by 

311 # ScalarFunction 

312 if bounds is None or len(bounds) == 0: 

313 new_bounds = (-np.inf, np.inf) 

314 else: 

315 new_bounds = old_bound_to_new(bounds) 

316 

317 # clip the initial guess to bounds, otherwise ScalarFunction doesn't work 

318 x = np.clip(x, new_bounds[0], new_bounds[1]) 

319 

320 # Set the parameters that SLSQP will need 

321 # meq, mieq: number of equality and inequality constraints 

322 meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args'])) 

323 for c in cons['eq']])) 

324 mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args'])) 

325 for c in cons['ineq']])) 

326 # m = The total number of constraints 

327 m = meq + mieq 

328 # la = The number of constraints, or 1 if there are no constraints 

329 la = array([1, m]).max() 

330 # n = The number of independent variables 

331 n = len(x) 

332 

333 # Define the workspaces for SLSQP 

334 n1 = n + 1 

335 mineq = m - meq + n1 + n1 

336 len_w = (3*n1+m)*(n1+1)+(n1-meq+1)*(mineq+2) + 2*mineq+(n1+mineq)*(n1-meq) \ 

337 + 2*meq + n1 + ((n+1)*n)//2 + 2*m + 3*n + 3*n1 + 1 

338 len_jw = mineq 

339 w = zeros(len_w) 

340 jw = zeros(len_jw) 

341 

342 # Decompose bounds into xl and xu 

343 if bounds is None or len(bounds) == 0: 

344 xl = np.empty(n, dtype=float) 

345 xu = np.empty(n, dtype=float) 

346 xl.fill(np.nan) 

347 xu.fill(np.nan) 

348 else: 

349 bnds = array(bounds, float) 

350 if bnds.shape[0] != n: 

351 raise IndexError('SLSQP Error: the length of bounds is not ' 

352 'compatible with that of x0.') 

353 

354 with np.errstate(invalid='ignore'): 

355 bnderr = bnds[:, 0] > bnds[:, 1] 

356 

357 if bnderr.any(): 

358 raise ValueError('SLSQP Error: lb > ub in bounds %s.' % 

359 ', '.join(str(b) for b in bnderr)) 

360 xl, xu = bnds[:, 0], bnds[:, 1] 

361 

362 # Mark infinite bounds with nans; the Fortran code understands this 

363 infbnd = ~isfinite(bnds) 

364 xl[infbnd[:, 0]] = np.nan 

365 xu[infbnd[:, 1]] = np.nan 

366 

367 # ScalarFunction provides function and gradient evaluation 

368 sf = _prepare_scalar_function(func, x, jac=jac, args=args, epsilon=eps, 

369 finite_diff_rel_step=finite_diff_rel_step, 

370 bounds=new_bounds) 

371 

372 # Initialize the iteration counter and the mode value 

373 mode = array(0, int) 

374 acc = array(acc, float) 

375 majiter = array(iter, int) 

376 majiter_prev = 0 

377 

378 # Initialize internal SLSQP state variables 

379 alpha = array(0, float) 

380 f0 = array(0, float) 

381 gs = array(0, float) 

382 h1 = array(0, float) 

383 h2 = array(0, float) 

384 h3 = array(0, float) 

385 h4 = array(0, float) 

386 t = array(0, float) 

387 t0 = array(0, float) 

388 tol = array(0, float) 

389 iexact = array(0, int) 

390 incons = array(0, int) 

391 ireset = array(0, int) 

392 itermx = array(0, int) 

393 line = array(0, int) 

394 n1 = array(0, int) 

395 n2 = array(0, int) 

396 n3 = array(0, int) 

397 

398 # Print the header if iprint >= 2 

399 if iprint >= 2: 

400 print("%5s %5s %16s %16s" % ("NIT", "FC", "OBJFUN", "GNORM")) 

401 

402 # mode is zero on entry, so call objective, constraints and gradients 

403 # there should be no func evaluations here because it's cached from 

404 # ScalarFunction 

405 fx = sf.fun(x) 

406 try: 

407 fx = float(np.asarray(fx)) 

408 except (TypeError, ValueError): 

409 raise ValueError("Objective function must return a scalar") 

410 g = append(sf.grad(x), 0.0) 

411 c = _eval_constraint(x, cons) 

412 a = _eval_con_normals(x, cons, la, n, m, meq, mieq) 

413 

414 while 1: 

415 # Call SLSQP 

416 slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw, 

417 alpha, f0, gs, h1, h2, h3, h4, t, t0, tol, 

418 iexact, incons, ireset, itermx, line, 

419 n1, n2, n3) 

420 

421 if mode == 1: # objective and constraint evaluation required 

422 fx = sf.fun(x) 

423 c = _eval_constraint(x, cons) 

424 

425 if mode == -1: # gradient evaluation required 

426 g = append(sf.grad(x), 0.0) 

427 a = _eval_con_normals(x, cons, la, n, m, meq, mieq) 

428 

429 if majiter > majiter_prev: 

430 # call callback if major iteration has incremented 

431 if callback is not None: 

432 callback(np.copy(x)) 

433 

434 # Print the status of the current iterate if iprint > 2 

435 if iprint >= 2: 

436 print("%5i %5i % 16.6E % 16.6E" % (majiter, sf.nfev, 

437 fx, linalg.norm(g))) 

438 

439 # If exit mode is not -1 or 1, slsqp has completed 

440 if abs(mode) != 1: 

441 break 

442 

443 majiter_prev = int(majiter) 

444 

445 # Optimization loop complete. Print status if requested 

446 if iprint >= 1: 

447 print(exit_modes[int(mode)] + " (Exit mode " + str(mode) + ')') 

448 print(" Current function value:", fx) 

449 print(" Iterations:", majiter) 

450 print(" Function evaluations:", sf.nfev) 

451 print(" Gradient evaluations:", sf.ngev) 

452 

453 return OptimizeResult(x=x, fun=fx, jac=g[:-1], nit=int(majiter), 

454 nfev=sf.nfev, njev=sf.ngev, status=int(mode), 

455 message=exit_modes[int(mode)], success=(mode == 0)) 

456 

457 

458def _eval_constraint(x, cons): 

459 # Compute constraints 

460 if cons['eq']: 

461 c_eq = concatenate([atleast_1d(con['fun'](x, *con['args'])) 

462 for con in cons['eq']]) 

463 else: 

464 c_eq = zeros(0) 

465 

466 if cons['ineq']: 

467 c_ieq = concatenate([atleast_1d(con['fun'](x, *con['args'])) 

468 for con in cons['ineq']]) 

469 else: 

470 c_ieq = zeros(0) 

471 

472 # Now combine c_eq and c_ieq into a single matrix 

473 c = concatenate((c_eq, c_ieq)) 

474 return c 

475 

476 

477def _eval_con_normals(x, cons, la, n, m, meq, mieq): 

478 # Compute the normals of the constraints 

479 if cons['eq']: 

480 a_eq = vstack([con['jac'](x, *con['args']) 

481 for con in cons['eq']]) 

482 else: # no equality constraint 

483 a_eq = zeros((meq, n)) 

484 

485 if cons['ineq']: 

486 a_ieq = vstack([con['jac'](x, *con['args']) 

487 for con in cons['ineq']]) 

488 else: # no inequality constraint 

489 a_ieq = zeros((mieq, n)) 

490 

491 # Now combine a_eq and a_ieq into a single a matrix 

492 if m == 0: # no constraints 

493 a = zeros((la, n)) 

494 else: 

495 a = vstack((a_eq, a_ieq)) 

496 a = concatenate((a, zeros([la, 1])), 1) 

497 

498 return a 

499 

500 

501if __name__ == '__main__': 

502 

503 # objective function 

504 def fun(x, r=[4, 2, 4, 2, 1]): 

505 """ Objective function """ 

506 return exp(x[0]) * (r[0] * x[0]**2 + r[1] * x[1]**2 + 

507 r[2] * x[0] * x[1] + r[3] * x[1] + 

508 r[4]) 

509 

510 # bounds 

511 bnds = array([[-inf]*2, [inf]*2]).T 

512 bnds[:, 0] = [0.1, 0.2] 

513 

514 # constraints 

515 def feqcon(x, b=1): 

516 """ Equality constraint """ 

517 return array([x[0]**2 + x[1] - b]) 

518 

519 def jeqcon(x, b=1): 

520 """ Jacobian of equality constraint """ 

521 return array([[2*x[0], 1]]) 

522 

523 def fieqcon(x, c=10): 

524 """ Inequality constraint """ 

525 return array([x[0] * x[1] + c]) 

526 

527 def jieqcon(x, c=10): 

528 """ Jacobian of inequality constraint """ 

529 return array([[1, 1]]) 

530 

531 # constraints dictionaries 

532 cons = ({'type': 'eq', 'fun': feqcon, 'jac': jeqcon, 'args': (1, )}, 

533 {'type': 'ineq', 'fun': fieqcon, 'jac': jieqcon, 'args': (10,)}) 

534 

535 # Bounds constraint problem 

536 print(' Bounds constraints '.center(72, '-')) 

537 print(' * fmin_slsqp') 

538 x, f = fmin_slsqp(fun, array([-1, 1]), bounds=bnds, disp=1, 

539 full_output=True)[:2] 

540 print(' * _minimize_slsqp') 

541 res = _minimize_slsqp(fun, array([-1, 1]), bounds=bnds, 

542 **{'disp': True}) 

543 

544 # Equality and inequality constraints problem 

545 print(' Equality and inequality constraints '.center(72, '-')) 

546 print(' * fmin_slsqp') 

547 x, f = fmin_slsqp(fun, array([-1, 1]), 

548 f_eqcons=feqcon, fprime_eqcons=jeqcon, 

549 f_ieqcons=fieqcon, fprime_ieqcons=jieqcon, 

550 disp=1, full_output=True)[:2] 

551 print(' * _minimize_slsqp') 

552 res = _minimize_slsqp(fun, array([-1, 1]), constraints=cons, 

553 **{'disp': True})