forked from MATPOWER/matpower
-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathqps_mosek.m
382 lines (370 loc) · 13 KB
/
qps_mosek.m
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
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
function [x, f, eflag, output, lambda] = qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt)
%QPS_MOSEK Quadratic Program Solver based on MOSEK.
% [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
% QPS_MOSEK(H, C, A, L, U, XMIN, XMAX, X0, OPT)
% A wrapper function providing a MATPOWER standardized interface for using
% MOSEKOPT 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)
%
% Inputs (all optional except H, C, A and L):
% H : matrix (possibly sparse) of quadratic cost coefficients
% C : vector of linear cost coefficients
% A, L, U : define the optional linear constraints. Default
% values for the elements of L and U are -Inf and Inf,
% respectively.
% XMIN, XMAX : optional lower and upper bounds on the
% X variables, defaults are -Inf and Inf, respectively.
% X0 : optional starting value of optimization vector X
% OPT : optional options structure with the following fields,
% all of which are also optional (default values shown in
% parentheses)
% verbose (0) - controls level of progress output displayed
% 0 = no progress output
% 1 = some progress output
% 2 = verbose progress output
% mosek_opt - options struct for MOSEK, value in verbose
% overrides these options
% PROBLEM : The inputs can alternatively be supplied in a single
% PROBLEM struct with fields corresponding to the input arguments
% described above: H, c, A, l, u, xmin, xmax, x0, opt
%
% Outputs:
% X : solution vector
% F : final objective function value
% EXITFLAG : exit flag
% 1 = success
% 0 = terminated at maximum number of iterations
% -1 = primal or dual infeasible
% < 0 = the negative of the MOSEK return code
% OUTPUT : output struct with the following fields:
% r - MOSEK return code
% res - MOSEK result struct
% LAMBDA : struct containing the Langrange and Kuhn-Tucker
% multipliers on the constraints, with fields:
% mu_l - lower (left-hand) limit on linear constraints
% mu_u - upper (right-hand) limit on linear constraints
% lower - lower bound on optimization variables
% upper - upper bound on optimization variables
%
% 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 A, L, U instead of
% A, B, Aeq, Beq.
%
% Calling syntax options:
% [x, f, exitflag, output, lambda] = ...
% qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt)
%
% x = qps_mosek(H, c, A, l, u)
% x = qps_mosek(H, c, A, l, u, xmin, xmax)
% x = qps_mosek(H, c, A, l, u, xmin, xmax, x0)
% x = qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt)
% x = qps_mosek(problem), where problem is a struct with fields:
% H, c, A, l, u, xmin, xmax, x0, opt
% all fields except 'c', 'A' and 'l' or 'u' are optional
% x = qps_mosek(...)
% [x, f] = qps_mosek(...)
% [x, f, exitflag] = qps_mosek(...)
% [x, f, exitflag, output] = qps_mosek(...)
% [x, f, exitflag, output, lambda] = qps_mosek(...)
%
% Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm)
% H = [ 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,1);
% A = [ 1 1 1 1;
% 0.17 0.11 0.10 0.18 ];
% l = [1; 0.10];
% u = [1; Inf];
% xmin = zeros(4,1);
% x0 = [1; 0; 0; 1];
% opt = struct('verbose', 2);
% [x, f, s, out, lambda] = qps_mosek(H, c, A, l, u, xmin, [], x0, opt);
%
% See also MOSEKOPT.
% MATPOWER
% Copyright (c) 2010-2016, Power Systems Engineering Research Center (PSERC)
% by Ray Zimmerman, PSERC Cornell
%
% This file is part of MATPOWER.
% Covered by the 3-clause BSD License (see LICENSE file for details).
% See http://www.pserc.cornell.edu/matpower/ for more info.
%% check for Optimization Toolbox
% if ~have_fcn('mosek')
% error('qps_mosek: requires MOSEK');
% end
%%----- input argument handling -----
%% gather inputs
if nargin == 1 && isstruct(H) %% problem struct
p = H;
else %% individual args
p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u);
if nargin > 5
p.xmin = xmin;
if nargin > 6
p.xmax = xmax;
if nargin > 7
p.x0 = x0;
if nargin > 8
p.opt = opt;
end
end
end
end
end
%% define nx, set default values for H and c
if ~isfield(p, 'H') || isempty(p.H) || ~any(any(p.H))
if (~isfield(p, 'A') || isempty(p.A)) && ...
(~isfield(p, 'xmin') || isempty(p.xmin)) && ...
(~isfield(p, 'xmax') || isempty(p.xmax))
error('qps_mosek: LP problem must include constraints or variable bounds');
else
if isfield(p, 'A') && ~isempty(p.A)
nx = size(p.A, 2);
elseif isfield(p, 'xmin') && ~isempty(p.xmin)
nx = length(p.xmin);
else % if isfield(p, 'xmax') && ~isempty(p.xmax)
nx = length(p.xmax);
end
end
p.H = sparse(nx, nx);
qp = 0;
else
nx = size(p.H, 1);
qp = 1;
end
if ~isfield(p, 'c') || isempty(p.c)
p.c = zeros(nx, 1);
end
if ~isfield(p, 'x0') || isempty(p.x0)
p.x0 = zeros(nx, 1);
end
%% default options
if ~isfield(p, 'opt')
p.opt = [];
end
if ~isempty(p.opt) && isfield(p.opt, 'verbose') && ~isempty(p.opt.verbose)
verbose = p.opt.verbose;
else
verbose = 0;
end
if ~isempty(p.opt) && isfield(p.opt, 'mosek_opt') && ~isempty(p.opt.mosek_opt)
mosek_opt = mosek_options(p.opt.mosek_opt);
else
mosek_opt = mosek_options;
end
%% set up problem struct for MOSEK
prob.c = p.c;
if qp
[prob.qosubi, prob.qosubj, prob.qoval] = find(tril(sparse(p.H)));
end
if isfield(p, 'A') && ~isempty(p.A)
prob.a = sparse(p.A);
nA = size(p.A, 1);
else
nA = 0;
end
if isfield(p, 'l') && ~isempty(p.A)
prob.blc = p.l;
end
if isfield(p, 'u') && ~isempty(p.A)
prob.buc = p.u;
end
if isfield(p, 'xmin') && ~isempty(p.xmin)
prob.blx = p.xmin;
end
if isfield(p, 'xmax') && ~isempty(p.xmax)
prob.bux = p.xmax;
end
%% A is not allowed to be empty
if ~isfield(prob, 'a') || isempty(prob.a)
unconstrained = 1;
prob.a = sparse(1, 1, 1, 1, nx);
prob.blc = -Inf;
prob.buc = Inf;
else
unconstrained = 0;
end
%%----- run optimization -----
if verbose
s = have_fcn('mosek', 'all');
if s.vnum < 7
alg_names = { %% version 6.x
'default', %% 0 : MSK_OPTIMIZER_FREE
'interior point', %% 1 : MSK_OPTIMIZER_INTPNT
'<conic>', %% 2 : MSK_OPTIMIZER_CONIC
'<qcone>', %% 3 : MSK_OPTIMIZER_QCONE
'primal simplex', %% 4 : MSK_OPTIMIZER_PRIMAL_SIMPLEX
'dual simplex', %% 5 : MSK_OPTIMIZER_DUAL_SIMPLEX
'primal dual simplex', %% 6 : MSK_OPTIMIZER_PRIMAL_DUAL_SIMPLEX
'automatic simplex', %% 7 : MSK_OPTIMIZER_FREE_SIMPLEX
'<mixed int>', %% 8 : MSK_OPTIMIZER_MIXED_INT
'<nonconvex>', %% 9 : MSK_OPTIMIZER_NONCONVEX
'concurrent' %% 10 : MSK_OPTIMIZER_CONCURRENT
};
elseif s.vnum < 8
alg_names = { %% version 7.x
'default', %% 0 : MSK_OPTIMIZER_FREE
'interior point', %% 1 : MSK_OPTIMIZER_INTPNT
'<conic>', %% 2 : MSK_OPTIMIZER_CONIC
'primal simplex', %% 3 : MSK_OPTIMIZER_PRIMAL_SIMPLEX
'dual simplex', %% 4 : MSK_OPTIMIZER_DUAL_SIMPLEX
'primal dual simplex', %% 5 : MSK_OPTIMIZER_PRIMAL_DUAL_SIMPLEX
'automatic simplex', %% 6 : MSK_OPTIMIZER_FREE_SIMPLEX
'network simplex', %% 7 : MSK_OPTIMIZER_NETWORK_PRIMAL_SIMPLEX
'<mixed int conic>', %% 8 : MSK_OPTIMIZER_MIXED_INT_CONIC
'<mixed int>', %% 9 : MSK_OPTIMIZER_MIXED_INT
'concurrent', %% 10 : MSK_OPTIMIZER_CONCURRENT
'<nonconvex>' %% 11 : MSK_OPTIMIZER_NONCONVEX
};
else
alg_names = { %% version 8.x
'<conic>', %% 0 : MSK_OPTIMIZER_CONIC
'dual simplex', %% 1 : MSK_OPTIMIZER_DUAL_SIMPLEX
'default', %% 2 : MSK_OPTIMIZER_FREE
'automatic simplex', %% 3 : MSK_OPTIMIZER_FREE_SIMPLEX
'interior point', %% 4 : MSK_OPTIMIZER_INTPNT
'<mixed int>', %% 5 : MSK_OPTIMIZER_MIXED_INT
'primal simplex' %% 6 : MSK_OPTIMIZER_PRIMAL_SIMPLEX
};
end
if qp
lpqp = 'QP';
else
lpqp = 'LP';
end
vn = have_fcn('mosek', 'vstr');
if isempty(vn)
vn = '<unknown>';
end
fprintf('MOSEK Version %s -- %s %s solver\n', ...
vn, alg_names{mosek_opt.MSK_IPAR_OPTIMIZER+1}, lpqp);
end
cmd = sprintf('minimize echo(%d)', verbose);
[r, res] = mosekopt(cmd, prob, mosek_opt);
%%----- repackage results -----
if isfield(res, 'sol')
if isfield(res.sol, 'bas')
sol = res.sol.bas;
else
sol = res.sol.itr;
end
x = sol.xx;
else
sol = [];
x = NaN(nx, 1);
end
%%----- process return codes -----
if isfield(res, 'symbcon')
sc = res.symbcon;
else
sc = mosek_symbcon;
end
eflag = -r;
msg = '';
switch (r)
case sc.MSK_RES_OK
if ~isempty(sol)
% if sol.solsta == sc.MSK_SOL_STA_OPTIMAL
if strcmp(sol.solsta, 'OPTIMAL')
msg = 'The solution is optimal.';
eflag = 1;
else
eflag = -1;
% if sol.prosta == sc.MSK_PRO_STA_PRIM_INFEAS
if strcmp(sol.prosta, 'PRIMAL_INFEASIBLE')
msg = 'The problem is primal infeasible.';
% elseif sol.prosta == sc.MSK_PRO_STA_DUAL_INFEAS
elseif strcmp(sol.prosta, 'DUAL_INFEASIBLE')
msg = 'The problem is dual infeasible.';
else
msg = sol.solsta;
end
end
end
case sc.MSK_RES_TRM_MAX_ITERATIONS
eflag = 0;
msg = 'The optimizer terminated at the maximum number of iterations.';
otherwise
if isfield(res, 'rmsg') && isfield(res, 'rcodestr')
msg = sprintf('%s : %s', res.rcodestr, res.rmsg);
else
msg = sprintf('MOSEK return code = %d', r);
end
end
if (verbose || r == sc.MSK_RES_ERR_LICENSE || ...
r == sc.MSK_RES_ERR_LICENSE_EXPIRED || ...
r == sc.MSK_RES_ERR_LICENSE_VERSION || ...
r == sc.MSK_RES_ERR_LICENSE_NO_SERVER_SUPPORT || ...
r == sc.MSK_RES_ERR_LICENSE_FEATURE || ...
r == sc.MSK_RES_ERR_LICENSE_INVALID_HOSTID || ...
r == sc.MSK_RES_ERR_LICENSE_SERVER_VERSION || ...
r == sc.MSK_RES_ERR_MISSING_LICENSE_FILE) ...
&& ~isempty(msg) %% always alert user of license problems
fprintf('%s\n', msg);
end
%%----- repackage results -----
if nargout > 1
if r == 0
f = p.c' * x;
if ~isempty(p.H)
f = 0.5 * x' * p.H * x + f;
end
else
f = [];
end
if nargout > 3
output.r = r;
output.res = res;
if nargout > 4
if ~isempty(sol)
if isfield(sol, 'slx')
lambda.lower = sol.slx;
else
lambda.lower = [];
end
if isfield(sol, 'sux')
lambda.upper = sol.sux;
else
lambda.upper = [];
end
if isfield(sol, 'slc')
lambda.mu_l = sol.slc;
else
lambda.mu_l = [];
end
if isfield(sol, 'suc')
lambda.mu_u = sol.suc;
else
lambda.mu_u = [];
end
else
if isfield(p, 'xmin') && ~isempty(p.xmin)
lambda.lower = NaN(nx, 1);
else
lambda.lower = [];
end
if isfield(p, 'xmax') && ~isempty(p.xmax)
lambda.upper = NaN(nx, 1);
else
lambda.upper = [];
end
lambda.mu_l = NaN(nA, 1);
lambda.mu_u = NaN(nA, 1);
end
if unconstrained
lambda.mu_l = [];
lambda.mu_u = [];
end
end
end
end