1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 """Quadratic Program Solver for PYPOWER.
18 """
19
20 import sys
21
22 from pypower.qps_pips import qps_pips
23 from pypower.qps_ipopt import qps_ipopt
24 from pypower.qps_cplex import qps_cplex
25 from pypower.qps_mosek import qps_mosek
26
27
28 -def qps_pypower(H, c=None, A=None, l=None, u=None, xmin=None, xmax=None,
29 x0=None, opt=None):
30 """Quadratic Program Solver for PYPOWER.
31
32 A common wrapper function for various QP solvers.
33 Solves the following QP (quadratic programming) problem::
34
35 min 1/2 x'*H*x + c'*x
36 x
37
38 subject to::
39
40 l <= A*x <= u (linear constraints)
41 xmin <= x <= xmax (variable bounds)
42
43 Inputs (all optional except C{H}, C{c}, C{A} and C{l}):
44 - C{H} : matrix (possibly sparse) of quadratic cost coefficients
45 - C{c} : vector of linear cost coefficients
46 - C{A, l, u} : define the optional linear constraints. Default
47 values for the elements of C{l} and C{u} are -Inf and Inf,
48 respectively.
49 - C{xmin}, C{xmax} : optional lower and upper bounds on the
50 C{x} variables, defaults are -Inf and Inf, respectively.
51 - C{x0} : optional starting value of optimization vector C{x}
52 - C{opt} : optional options structure with the following fields,
53 all of which are also optional (default values shown in parentheses)
54 - C{alg} (0) - determines which solver to use
55 - 0 = automatic, first available of BPMPD_MEX, CPLEX, MIPS
56 - 100 = BPMPD_MEX
57 - 200 = MIPS, MATLAB Interior Point Solver
58 pure MATLAB implementation of a primal-dual
59 interior point method
60 - 250 = MIPS-sc, a step controlled variant of MIPS
61 - 300 = Optimization Toolbox, QUADPROG or LINPROG
62 - 400 = IPOPT
63 - 500 = CPLEX
64 - 600 = MOSEK
65 - C{verbose} (0) - controls level of progress output displayed
66 - 0 = no progress output
67 - 1 = some progress output
68 - 2 = verbose progress output
69 - C{max_it} (0) - maximum number of iterations allowed
70 - 0 = use algorithm default
71 - C{bp_opt} - options vector for BP
72 - C{cplex_opt} - options dict for CPLEX
73 - C{ipopt_opt} - options dict for IPOPT
74 - C{pips_opt} - options dict for L{qps_pips}
75 - C{mosek_opt} - options dict for MOSEK
76 - C{ot_opt} - options dict for QUADPROG/LINPROG
77 - C{problem} : The inputs can alternatively be supplied in a single
78 C{problem} dict with fields corresponding to the input arguments
79 described above: C{H, c, A, l, u, xmin, xmax, x0, opt}
80
81 Outputs:
82 - C{x} : solution vector
83 - C{f} : final objective function value
84 - C{exitflag} : exit flag
85 - 1 = converged
86 - 0 or negative values = algorithm specific failure codes
87 - C{output} : output struct with the following fields:
88 - C{alg} - algorithm code of solver used
89 - (others) - algorithm specific fields
90 - C{lmbda} : dict containing the Langrange and Kuhn-Tucker
91 multipliers on the constraints, with fields:
92 - C{mu_l} - lower (left-hand) limit on linear constraints
93 - C{mu_u} - upper (right-hand) limit on linear constraints
94 - C{lower} - lower bound on optimization variables
95 - C{upper} - upper bound on optimization variables
96
97
98 Example from U{http://www.uc.edu/sashtml/iml/chap8/sect12.htm}:
99
100 >>> from numpy import array, zeros, Inf
101 >>> from scipy.sparse import csr_matrix
102 >>> H = csr_matrix(array([[1003.1, 4.3, 6.3, 5.9],
103 ... [4.3, 2.2, 2.1, 3.9],
104 ... [6.3, 2.1, 3.5, 4.8],
105 ... [5.9, 3.9, 4.8, 10 ]]))
106 >>> c = zeros(4)
107 >>> A = csr_matrix(array([[1, 1, 1, 1 ],
108 ... [0.17, 0.11, 0.10, 0.18]]))
109 >>> l = array([1, 0.10])
110 >>> u = array([1, Inf])
111 >>> xmin = zeros(4)
112 >>> xmax = None
113 >>> x0 = array([1, 0, 0, 1])
114 >>> solution = qps_pips(H, c, A, l, u, xmin, xmax, x0)
115 >>> round(solution["f"], 11) == 1.09666678128
116 True
117 >>> solution["converged"]
118 True
119 >>> solution["output"]["iterations"]
120 10
121
122 @author: Ray Zimmerman (PSERC Cornell)
123 @author: Richard Lincoln
124 """
125
126
127 if isinstance(H, dict):
128 p = H
129 if 'opt' in p: opt = p['opt']
130 if 'x0' in p: x0 = p['x0']
131 if 'xmax' in p: xmax = p['xmax']
132 if 'xmin' in p: xmin = p['xmin']
133 if 'u' in p: u = p['u']
134 if 'l' in p: l = p['l']
135 if 'A' in p: A = p['A']
136 if 'c' in p: c = p['c']
137 if 'H' in p: H = p['H']
138 else:
139
140 assert c is not None
141
142
143
144 if opt is None:
145 opt = {}
146
147
148
149
150
151
152
153
154 if 'alg' in opt:
155 alg = opt['alg']
156 else:
157 alg = 0
158
159 if 'verbose' in opt:
160 verbose = opt['verbose']
161 else:
162 verbose = 0
163
164 if alg == 0:
165 try:
166 import pymosek
167 alg = 600
168 except ImportError:
169 try:
170 import cplex
171 alg = 500
172 except ImportError:
173 alg = 200
174
175
176
177 if alg == 200 or alg == 250:
178
179 if 'pips_opt' in opt:
180 pips_opt = opt['pips_opt']
181 else:
182 pips_opt = {}
183
184 if 'max_it' in opt:
185 pips_opt['max_it'] = opt['max_it']
186
187 if alg == 200:
188 pips_opt['step_control'] = False
189 else:
190 pips_opt['step_control'] = True
191
192 pips_opt['verbose'] = verbose
193
194
195 x, f, eflag, output, lmbda = \
196 qps_pips(H, c, A, l, u, xmin, xmax, x0, pips_opt)
197 elif alg == 400:
198 x, f, eflag, output, lmbda = \
199 qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt)
200 elif alg == 500:
201 x, f, eflag, output, lmbda = \
202 qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt)
203 elif alg == 600:
204 x, f, eflag, output, lmbda = \
205 qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt)
206 else:
207 sys.stderr.write('qps_pypower: %d is not a valid algorithm code\n', alg)
208
209 if 'alg' not in output:
210 output['alg'] = alg
211
212 return x, f, eflag, output, lmbda
213
214
215 if __name__ == "__main__":
216 import doctest
217 doctest.testmod()
218