-
Notifications
You must be signed in to change notification settings - Fork 465
/
qps_pips.py
179 lines (153 loc) · 7.07 KB
/
qps_pips.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
# Copyright (c) 1996-2015 PSERC. All rights reserved.
# Use of this source code is governed by a BSD-style
# license that can be found in the LICENSE file.
"""Uses the Python Interior Point Solver (PIPS) to solve QP (quadratic
programming) problems.
"""
from numpy import Inf, ones, zeros, dot
from scipy.sparse import csr_matrix as sparse
from pandapower.pypower.pips import pips
def qps_pips(H, c, A, l, u, xmin=None, xmax=None, x0=None, opt=None):
"""Uses the Python Interior Point Solver (PIPS) to solve the following
QP (quadratic programming) problem::
min 1/2 x'*H*x + C'*x
x
subject to::
l <= A*x <= u (linear constraints)
xmin <= x <= xmax (variable bounds)
Note the calling syntax is almost identical to that of QUADPROG from
MathWorks' Optimization Toolbox. The main difference is that the linear
constraints are specified with C{A}, C{L}, C{U} instead of C{A}, C{B},
C{Aeq}, C{Beq}.
Example from U{http://www.uc.edu/sashtml/iml/chap8/sect12.htm}:
>>> from numpy import array, zeros, Inf
>>> from scipy.sparse import csr_matrix
>>> H = csr_matrix(array([[1003.1, 4.3, 6.3, 5.9],
... [4.3, 2.2, 2.1, 3.9],
... [6.3, 2.1, 3.5, 4.8],
... [5.9, 3.9, 4.8, 10 ]]))
>>> c = zeros(4)
>>> A = csr_matrix(array([[1, 1, 1, 1 ],
... [0.17, 0.11, 0.10, 0.18]]))
>>> l = array([1, 0.10])
>>> u = array([1, Inf])
>>> xmin = zeros(4)
>>> xmax = None
>>> x0 = array([1, 0, 0, 1])
>>> solution = qps_pips(H, c, A, l, u, xmin, xmax, x0)
>>> round(solution["f"], 11) == 1.09666678128
True
>>> solution["converged"]
True
>>> solution["output"]["iterations"]
10
All parameters are optional except C{H}, C{c}, C{A} and C{l} or C{u}.
@param H: Quadratic cost coefficients.
@type H: csr_matrix
@param c: vector of linear cost coefficients
@type c: array
@param A: Optional linear constraints.
@type A: csr_matrix
@param l: Optional linear constraints. Default values are M{-Inf}.
@type l: array
@param u: Optional linear constraints. Default values are M{Inf}.
@type u: array
@param xmin: Optional lower bounds on the M{x} variables, defaults are
M{-Inf}.
@type xmin: array
@param xmax: Optional upper bounds on the M{x} variables, defaults are
M{Inf}.
@type xmax: array
@param x0: Starting value of optimization vector M{x}.
@type x0: array
@param opt: optional options dictionary with the following keys, all of
which are also optional (default values shown in parentheses)
- C{verbose} (False) - Controls level of progress output
displayed
- C{feastol} (1e-6) - termination tolerance for feasibility
condition
- C{gradtol} (1e-6) - termination tolerance for gradient
condition
- C{comptol} (1e-6) - termination tolerance for
complementarity condition
- C{costtol} (1e-6) - termination tolerance for cost
condition
- C{max_it} (150) - maximum number of iterations
- C{step_control} (False) - set to True to enable step-size
control
- C{max_red} (20) - maximum number of step-size reductions if
step-control is on
- C{cost_mult} (1.0) - cost multiplier used to scale the
objective function for improved conditioning. Note: The
same value must also be passed to the Hessian evaluation
function so that it can appropriately scale the objective
function term in the Hessian of the Lagrangian.
@type opt: dict
@rtype: dict
@return: The solution dictionary has the following keys:
- C{x} - solution vector
- C{f} - final objective function value
- C{converged} - exit status
- True = first order optimality conditions satisfied
- False = maximum number of iterations reached
- None = numerically failed
- C{output} - output dictionary with keys:
- C{iterations} - number of iterations performed
- C{hist} - dictionary of arrays with trajectories of the
following: feascond, gradcond, coppcond, costcond, gamma,
stepsize, obj, alphap, alphad
- C{message} - exit message
- C{lmbda} - dictionary containing the Langrange and Kuhn-Tucker
multipliers on the constraints, with keys:
- C{eqnonlin} - nonlinear equality constraints
- C{ineqnonlin} - nonlinear inequality constraints
- C{mu_l} - lower (left-hand) limit on linear constraints
- C{mu_u} - upper (right-hand) limit on linear constraints
- C{lower} - lower bound on optimization variables
- C{upper} - upper bound on optimization variables
@see: L{pips}
@author: Ray Zimmerman (PSERC Cornell)
"""
if isinstance(H, dict):
p = H
else:
p = {'H': H, 'c': c, 'A': A, 'l': l, 'u': u}
if xmin is not None: p['xmin'] = xmin
if xmax is not None: p['xmax'] = xmax
if x0 is not None: p['x0'] = x0
if opt is not None: p['opt'] = opt
if 'H' not in p or p['H'] == None:#p['H'].nnz == 0:
if p['A'] is None or p['A'].nnz == 0 and \
'xmin' not in p and \
'xmax' not in p:
# 'xmin' not in p or len(p['xmin']) == 0 and \
# 'xmax' not in p or len(p['xmax']) == 0:
print('qps_pips: LP problem must include constraints or variable bounds')
return
else:
if p['A'] is not None and p['A'].nnz >= 0:
nx = p['A'].shape[1]
elif 'xmin' in p and len(p['xmin']) > 0:
nx = p['xmin'].shape[0]
elif 'xmax' in p and len(p['xmax']) > 0:
nx = p['xmax'].shape[0]
p['H'] = sparse((nx, nx))
else:
nx = p['H'].shape[0]
p['xmin'] = -Inf * ones(nx) if 'xmin' not in p else p['xmin']
p['xmax'] = Inf * ones(nx) if 'xmax' not in p else p['xmax']
p['c'] = zeros(nx) if p['c'] is None else p['c']
p['x0'] = zeros(nx) if 'x0' not in p else p['x0']
def qp_f(x, return_hessian=False):
f = 0.5 * dot(x * p['H'], x) + dot(p['c'], x)
df = p['H'] * x + p['c']
if not return_hessian:
return f, df
d2f = p['H']
return f, df, d2f
p['f_fcn'] = qp_f
sol = pips(p)
return sol["x"], sol["f"], sol["eflag"], sol["output"], sol["lmbda"]
if __name__ == "__main__":
import doctest
doctest.testmod()