From 719457270c79d42f763da3a629dcb33d1b7a7811 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilson=20Gonz=C3=A1lez=20Vanegas?= Date: Fri, 31 Jul 2026 17:55:20 -0500 Subject: [PATCH 1/4] Add wrapper function for solving MIQCQP problems --- lib/miqcqps_master.m | 288 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 288 insertions(+) create mode 100644 lib/miqcqps_master.m diff --git a/lib/miqcqps_master.m b/lib/miqcqps_master.m new file mode 100644 index 0000000..4969c66 --- /dev/null +++ b/lib/miqcqps_master.m @@ -0,0 +1,288 @@ +function [x, f, eflag, output, lambda] = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt) +% miqcqps_master - Mixed Integer Quadratically Constrained Quadratic Program Solver wrapper function. +% :: +% +% [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... +% MIQCQPS_MASTER(H, C, Q, B, LQ, UQ, A, L, U, XMIN, XMAX, X0, VTYPE, OPT) +% [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQCQPS_MASTER(PROBLEM) +% A common wrapper function for various MIQCQP solvers. +% Solves the following MIQCQP (mixed integer quadratically constrained quadratic +% programming) problem: +% +% min 1/2 X'*H*X + C'*X +% X +% +% subject to +% +% LQ(i) <= 1/2 X'*Q{i}*X + B(i,:)*X <= UQ(i), i = 1,2,...,nq +% (quadratic constraints) +% L <= A*X <= U (linear constraints) +% XMIN <= X <= XMAX (variable bounds) +% X(i) is integer, for i in I (integer variable constraints) +% X(b) is binary, for b in B (binary variable constraints) +% +% Inputs (all optional except H, C, Q, B, LQ, and UQ): +% H : matrix (possibly sparse) of quadratic cost coefficients +% C : vector of linear cost coefficients +% Q : nq x 1 cell array of sparse quadratic matrices for quadratic constraints +% B : matrix (possibly sparse) of linear term of quadratic constraints +% LQ, UQ: define the lower an upper bounds on the quadratic constraints +% 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 +% VTYPE : character string of length NX (number of elements in X), +% or 1 (value applies to all variables in x), +% allowed values are 'C' (continuous), 'B' (binary), or +% 'I' (integer), 'S' (semi-continuous), or 'N' (semi-integer). +% OPT : optional options structure with the following fields, +% all of which are also optional (default values shown in +% parentheses) +% alg ('DEFAULT') : determines which solver to use, can be either +% a string (new-style) or a numerical alg code (old-style) +% 'DEFAULT' : equals 'GUROBI' in current implementation +% 'GUROBI' : Gurobi, requires Gurobi solver +% https://www.gurobi.com +% verbose (0) - controls level of progress output displayed +% 0 = no progress output +% 1 = some progress output +% 2 = verbose progress output +% fix_integer (0) - fix integer variables at value in x0, if true +% relax_integer (0) - relax integer constraints, if true +% grb_opt - options struct for GUROBI +% PROBLEM : The inputs can alternatively be supplied in a single +% PROBLEM struct with fields corresponding to the input arguments +% described above: H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt +% +% Outputs: +% X : solution vector +% F : final objective function value +% EXITFLAG : exit flag +% 1 = converged +% 0 or negative values = solver specific failure codes +% OUTPUT : output struct with the following fields: +% alg - algorithm code of solver used +% (others) - algorithm specific fields +% 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 +% mu_lq - lower (left-hand) limit on quadratic constraints +% mu_uq - upper (right-hand) limit on quadratic constraints +% lower - lower bound on optimization variables +% upper - upper bound on optimization variables +% +% Calling syntax options: +% [x, f, exitflag, output, lambda] = ... +% miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt) +% +% x = miqcqps_master(H, c, Q, B, lq, uq) +% x = miqcqps_master(H, c, Q, B, lq, uq, A, l, u) +% x = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax) +% x = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0) +% x = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype) +% x = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt) +% x = miqcqps_master(problem), where problem is a struct with fields: +% H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt +% all fields except 'c', are optional, and problem with +% linear costs must include constraints +% x = miqcqps_master(...) +% [x, f] = miqcqps_master(...) +% [x, f, exitflag] = miqcqps_master(...) +% [x, f, exitflag, output] = miqcqps_master(...) +% [x, f, exitflag, output, lambda] = miqcqps_master(...) +% +% Example: (problem from OPTI Toolbox, see: +% https://jonathancurrie.github.io/OPTI/examples/problem-types/miqcqp/) +% +% H = [1 -1; -1 2]; +% c = [-2 -6]'; +% Q = {}; B = []; lq = []; uq = []; +% A = [1 1; -1 2; 2 1]; +% u = [2 2 3]'; +% l = []; +% xmin = zeros(2,1); +% xmax = Inf(2,1); +% x0 = []; +% vtype = 'IC'; +% opt = struct('verbose', 2); +% [x, f, s, out, lam] = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt); + +% MP-Opt-Model +% Copyright (c) 2019-2026, Power Systems Engineering Research Center (PSERC) +% by Wilson Gonzalez Vanegas, Universidad Nacional de Colombia Sede Manizales +% and Ray Zimmerman, PSERC Cornell +% +% This file is part of MP-Opt-Model. +% Covered by the 3-clause BSD License (see LICENSE file for details). +% See https://github.com/MATPOWER/mp-opt-model for more info. + +%%----- input argument handling ----- +%% gather inputs +if nargin == 1 && isstruct(H) %% problem struct + p = H; + if isfield(p, 'opt'), opt = p.opt; else, opt = []; end + if isfield(p, 'vtype'), vtype = p.vtype;else, vtype = []; end + if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end + if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end + if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end + if isfield(p, 'u'), u = p.u; else, u = []; end + if isfield(p, 'l'), l = p.l; else, l = []; end + if isfield(p, 'A'), A = p.A; else, A = []; end + if isfield(p, 'uq'), uq = p.uq; else, uq = []; end + if isfield(p, 'lq'), lq = p.lq; else, lq = []; end + if isfield(p, 'B'), B = p.B; else, B = []; end + if isfield(p, 'Q'), Q = p.Q; else, Q = {}; end + if isfield(p, 'c'), c = p.c; else, c = []; end + if isfield(p, 'H'), H = p.H; else, H = []; end +else %% individual args + if nargin < 14 + opt = []; + if nargin < 13 + vtype = []; + if nargin < 12 + x0 = []; + if nargin < 11 + xmax = []; + if nargin < 10 + xmin = []; + if nargin < 7 + A = []; + l = []; + u = []; + end + end + end + end + end + end +end + +%% default options +if ~isempty(opt) && isfield(opt, 'alg') && ~isempty(opt.alg) + alg = opt.alg; + % convert integer codes to string values + if ~ischar(alg) + switch alg + case 700 + alg = 'GUROBI'; + otherwise + error('miqcqps_master: %d is not a valid algorithm code. the only available solver with MIQCQP interface is Gurobi (700)', alg); + end + end +else + alg = 'DEFAULT'; +end +if strcmp(alg, 'DEFAULT') + if have_feature('gurobi') + alg = 'GUROBI'; + end +end + +%% handle relax_integer and fix_integer options +done = false; +fix_integer_presolve = true; +if ~isempty(vtype) && (isfield(opt, 'relax_integer') && opt.relax_integer || ... + isfield(opt, 'fix_integer') && opt.fix_integer) + nx = length(x0); + if length(vtype) == 1 %% expand if necessary + vtype = char(vtype * ones(1, nx)); + end + j = (vtype == 'B' | vtype == 'I')'; + if isfield(opt, 'fix_integer') && opt.fix_integer + x0(j) = round(x0(j)); + if fix_integer_presolve + x = x0; + % 1) Reorganize quadratic constraints + if isempty(Q) + QQ = []; + BB = B; + Bxj = 0; + else + [nq, ncolQ] = size(Q); %% Check dimension of quadratic terms of quadratic constraints + if ~iscell(Q) || ncolQ ~= 1 + error('miqcqps_master: Q must be column vector cell array.') + end + sizeQi = cell2mat(cellfun(@(x)(size(x)), Q, 'UniformOutput', false)); + if any(sizeQi(:) - sizeQi(1)) + error('miqcqps_master: All matrices Q{i}, i=1,...,%d must be square of the same size.', nq) + end + + QQ = cellfun(@(w)(w(~j,~j)), Q, 'UniformOutput', false); + BB = cell2mat(cellfun(@(w)( x(j)' * (1/2*(w(~j,j)'+w(j,~j))) ), Q, 'UniformOutput', false)); + Bxj = cell2mat(cellfun(@(w)( 1/2 * x(j)' * w(j,j) * x(j) ), Q, 'UniformOutput', false)); + end + % 2) Reorganize linear constraints + Axj = A(:,j) * x(j); + + % 3) Reorganize objective function + cc = c(~j); + if isempty(H) + HH = []; + else + HH = H(~j, ~j); + cc = cc + 1/2*(H(~j,j)+H(j,~j)') * x(j); + end + [x(~j), f, eflag, output, lambda] = ... + qcqps_master(HH, cc, QQ, BB, lq-Bxj, uq-Bxj, A(:, ~j), l-Axj, u-Axj, ... + xmin(~j), xmax(~j), x(~j), opt); + f = f + c(j)' * x(j); + if ~isempty(H) + f = f + 1/2 * x(j)' * H(j,j) * x(j); + end + mu_lower = zeros(size(x)); + mu_upper = zeros(size(x)); + mu_lower(~j) = lambda.lower; + mu_upper(~j) = lambda.upper; + lambda.lower = mu_lower; + lambda.upper = mu_upper; + done = true; + else + %% fix integer variables + if ~isempty(j) + %% expand if necessary + if length(xmin) == 1 + xmin = xmin * ones(nx, 1); + elseif isempty(xmin) + xmin = -Inf(nx, 1); + end + if length(xmax) == 1 + xmax = xmax * ones(nx, 1); + elseif isempty(xmax) + xmax = Inf(nx, 1); + end + xmin(j) = x0(j); + xmax(j) = x0(j); + end + end + end + if ~done + vtype(j) = 'C'; + end +end + +miqp = 1; +if any(cell2mat(cellfun(@(x)(any(x(:))), Q, 'UniformOutput', false))) + miqp = 0; +end + +%%----- call the appropriate solver ----- +if ~done + if miqp && ~strcmp(alg, 'GUROBI') + [x, f, eflag, output, lambda] = ... + miqps_master(H, c, A, l, u, xmin, xmax, x0, vtype, opt); + else + switch alg + case 'GUROBI' + [x, f, eflag, output, lambda] = ... + miqcqps_gurobi(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt); + otherwise + error('miqcqps_master: currently, the only available solver with MIQCQP interface is Gurobi') + end + end +end +if ~isfield(output, 'alg') || isempty(output.alg) + output.alg = alg; +end From d796dc61e333e268f61370c4e898fea4568f2f1c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilson=20Gonz=C3=A1lez=20Vanegas?= Date: Fri, 31 Jul 2026 17:56:19 -0500 Subject: [PATCH 2/4] Add MIQCQP solver based on Gurobi. --- lib/miqcqps_gurobi.m | 515 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 515 insertions(+) create mode 100644 lib/miqcqps_gurobi.m diff --git a/lib/miqcqps_gurobi.m b/lib/miqcqps_gurobi.m new file mode 100644 index 0000000..d29b602 --- /dev/null +++ b/lib/miqcqps_gurobi.m @@ -0,0 +1,515 @@ +function [x, f, eflag, output, lambda] = miqcqps_gurobi(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt) +% miqcqps_gurobi - Mixed Integer Quadratically Constrained Quadratic Program Solver based on GUROBI. +% :: +% +% [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... +% MIQCQPS_GUROBI(H, C, Q, B, LQ, UQ, A, L, U, XMIN, XMAX, X0, OPT) +% [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQCQPS_GUROBI(PROBLEM) +% A wrapper function providing a standardized interface for using +% GUROBI to solve the following (possibly non-convex) MIQCQP (Mixed Integer +% quadratically constrained quadratic programming) problem: +% +% min 1/2 X'*H*X + C'*X +% X +% +% subject to +% +% LQ(i) <= 1/2 X'*Q{i}*X + B(i,:)*X <= UQ(i), i = 1,2,...,nq +% (quadratic constraints) +% L <= A*X <= U (linear constraints) +% XMIN <= X <= XMAX (variable bounds) +% X(i) is integer, for i in I (integer variable constraints) +% X(b) is binary, for b in B (binary variable constraints) +% +% Inputs (all optional except H, C, Q, B, LQ, and UQ): +% H : matrix (possibly sparse) of quadratic cost coefficients +% C : vector of linear cost coefficients +% Q : nq x 1 cell array of sparse quadratic matrices for quadratic constraints +% B : matrix (possibly sparse) of linear term of quadratic constraints +% LQ, UQ: define the lower an upper bounds on the quadratic constraints +% 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 +% VTYPE : character string of length NX (number of elements in X), +% or 1 (value applies to all variables in x), +% allowed values are 'C' (continuous), 'B' (binary), or +% 'I' (integer), 'S' (semi-continuous), or 'N' (semi-integer). +% 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 +% 3 = even more verbose progress output +% grb_opt - options struct for GUROBI, 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, Q, B, lq, uq, A, l, u, xmin, xmax, x0, opt +% +% Outputs: +% X : solution vector +% F : final objective function value +% EXITFLAG : exit flag +% 1 = converged +% 0 or negative values = solver specific failure codes +% OUTPUT : output struct with the following fields: +% alg - algorithm code of solver used +% (others) - algorithm specific fields +% 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 +% mu_lq - lower (left-hand) limit on quadratic constraints +% mu_uq - upper (right-hand) limit on quadratic constraints +% lower - lower bound on optimization variables +% upper - upper bound on optimization variables +% +% Calling syntax options: +% [x, f, exitflag, output, lambda] = ... +% miqcqps_gurobi(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt) +% +% x = miqcqps_gurobi(H, c, Q, B, lq, uq) +% x = miqcqps_gurobi(H, c, Q, B, lq, uq, A, l, u) +% x = miqcqps_gurobi(H, c, Q, B, lq, uq, A, l, u, xmin, xmax) +% x = miqcqps_gurobi(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0) +% x = miqcqps_gurobi(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype) +% x = miqcqps_gurobi(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt) +% x = miqcqps_gurobi(problem), where problem is a struct with fields: +% H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt +% all fields except 'c', are optional, and problem with +% linear costs must include constraints +% x = miqcqps_gurobi(...) +% [x, f] = miqcqps_gurobi(...) +% [x, f, exitflag] = miqcqps_gurobi(...) +% [x, f, exitflag, output] = miqcqps_gurobi(...) +% [x, f, exitflag, output, lambda] = miqcqps_gurobi(...) +% +% Example: (problem from OPTI Toolbox, see: +% https://jonathancurrie.github.io/OPTI/examples/problem-types/miqcqp/) +% +% H = [1 -1; -1 2]; +% c = [-2 -6]'; +% Q = {}; B = []; lq = []; uq = []; +% A = [1 1; -1 2; 2 1]; +% u = [2 2 3]'; +% l = []; +% xmin = zeros(2,1); +% xmax = Inf(2,1); +% x0 = []; +% opt = struct('verbose', 2); +% vtype = 'IC'; +% [x, f, s, out, lam] = miqcqps_gurobi(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt); +% +% See also qcqps_master, gurobi_options, gurobi. + +% MP-Opt-Model +% Copyright (c) 2019-2026, Power Systems Engineering Research Center (PSERC) +% by Wilson Gonzalez Vanegas, Universidad Nacional de Colombia Sede Manizales +% and Ray Zimmerman, PSERC Cornell +% +% This file is part of MP-Opt-Model. +% Covered by the 3-clause BSD License (see LICENSE file for details). +% See https://github.com/MATPOWER/mp-opt-model for more info. + +%%----- input argument handling ----- +%% gather inputs +if nargin == 1 && isstruct(H) %% problem struct + p = H; + if isfield(p, 'opt'), opt = p.opt; else, opt = []; end + if isfield(p, 'vtype'), vtype = p.vtype;else, vtype = []; end + if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end + if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end + if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end + if isfield(p, 'u'), u = p.u; else, u = []; end + if isfield(p, 'l'), l = p.l; else, l = []; end + if isfield(p, 'A'), A = p.A; else, A = []; end + if isfield(p, 'uq'), uq = p.uq; else, uq = []; end + if isfield(p, 'lq'), lq = p.lq; else, lq = []; end + if isfield(p, 'B'), B = p.B; else, B = []; end + if isfield(p, 'Q'), Q = p.Q; else, Q = {}; end + if isfield(p, 'c'), c = p.c; else, c = []; end + if isfield(p, 'H'), H = p.H; else, H = []; end +else %% individual args + if nargin < 14 + opt = []; + if nargin < 13 + vtype = []; + if nargin < 12 + x0 = []; + if nargin < 11 + xmax = []; + if nargin < 10 + xmin = []; + if nargin < 7 + A = []; + l = []; + u = []; + end + end + end + end + end + end +end + +%% define nx, nq, nlin, and set default values for missing optional inputs +if ~isempty(Q) + [nq, ncolQ] = size(Q); + if ~iscell(Q) || ncolQ ~= 1 + error('miqcqps_gurobi: Q must be column vector cell array.') + end + sizeQi = cell2mat(cellfun(@(x)(size(x)), Q, 'UniformOutput', false)); + if any(sizeQi(:) - sizeQi(1)) + error('miqcqps_gurobi: All matrices Q{i}, i=1,...,%d must be square of the same size.', nq) + end +elseif ~isempty(B) + nq = size(B, 1); +else + nq = 0; +end + +if ~isempty(H) + [nrowH, ncolH] = size(H); + if nrowH ~= ncolH + error('miqcqps_gurobi: H must be a square matrix.') + end + nx = nrowH; +else + if ~isempty(c) + nx = length(c); + else + if nq + nx = size(Q{1},2); + else + if ~isempty(B) + nx = size(B,2); + else + if ~isempty(A) + nx = size(A,2); + else + error('miqcqps_gurobi: inputs arguments H, c, Q, B, and A can not be all empty.') + end + end + end + end +end + +if isempty(H) && isempty(c) && isempty(B) && isempty(A) + error('miqcqps_gurobi: Problem is incomplete: H, c, B and A can not be all empty.') +end +if isempty(H) + H = sparse(nx,nx); +end +if isempty(c) + c = sparse(nx,1); +elseif length(c) ~= nx + error('miqcqps_gurobi: Dimension of c (%d) must be iqual to the number of variables (%d).', length(c), nx) +end + +if nq + if isempty(B) + B = sparse(nq,nx); + else + [nrowB, ncolB] = size(B); + if nrowB ~= nq || ncolB ~= nx + error('miqcqps_gurobi: Dimension of B (%dx%d) should be number of quad constraints times number of variables (%dx%d).', nrowB, ncolB, nq, nx); + end + end +else + Q = {}; + if ~isempty(B) + [~, ncolB] = size(B); + if ncolB ~= nx + error('miqcqps_gurobi: The number of columns of matrix B (%d) must be equal to the number of variables (%d).', ncolB, nx); + end + else + B = sparse(nq,nx); + if ~isempty(lq) || ~isempty(uq) + error('miqcqps_gurobi: No quadratic constraints were found. Inputs lq and uq should be empty.') + end + end +end + +if ~isempty(A) + [nlin, ncolA] = size(A); + if ncolA ~= nx + error('miqcqps_gurobi: The number of columns of matrix A (%d) must be equal to the number of variables (%d).', ncolA, nx); + end +else + nlin = 0; + A = sparse(nlin,nx); +end + +if isempty(uq) %% By default, quadratic inequalities are ... + uq = Inf(nq, 1); %% ... unbounded above and ... +elseif length(uq) ~= size(B,1) + error('miqcqps_gurobi: Dimension mismatch between uq, Q, and B.') +end +if isempty(lq) + lq = -Inf(nq, 1); %% ... unbounded below. +elseif length(lq) ~= size(B,1) + error('miqcqps_gurobi: Dimension mismatch between lq, Q, and B.') +end +if isempty(u) %% By default, linear inequalities are ... + u = Inf(nlin, 1); %% ... unbounded above and ... +elseif length(u) ~= nlin + error('miqcqps_gurobi: Dimension of u (%d) must be iqual to the number of linear constraints (%d).', length(u), nlin') +end +if isempty(l) + l = -Inf(nlin, 1); %% ... unbounded below. +elseif length(l) ~= nlin + error('miqcqps_gurobi: Dimension of l (%d) must be iqual to the number of linear constraints (%d).', length(l), nlin') +end +if isempty(xmin) %% By default, optimization variables are ... + xmin = -Inf(nx, 1); %% ... unbounded below and ... +elseif length(xmin) ~= nx + error('miqcqps_gurobi: Dimension of xmin (%d) must be iqual to the number of variables (%d).', length(xmin), nx') +end +if isempty(xmax) + xmax = Inf(nx, 1); %% ... unbounded above. +elseif length(xmax) ~= nx + error('miqcqps_gurobi: Dimension of xmax (%d) must be iqual to the number of variables (%d).', length(xmax), nx') +end +if isempty(x0) + x0 = zeros(nx, 1); +elseif length(x0) ~= nx + error('miqcqps_gurobi: Dimension of x0 (%d) must be iqual to the number of variables (%d).', length(x0), nx') +end + +%% default options +if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) + verbose = opt.verbose; +else + verbose = 0; +end + +%% set up options struct for Gurobi +if ~isempty(opt) && isfield(opt, 'grb_opt') && ~isempty(opt.grb_opt) + g_opt = gurobi_options(opt.grb_opt); +else + g_opt = gurobi_options; +end +if verbose > 1 + g_opt.LogToConsole = 1; + g_opt.OutputFlag = 1; + if verbose > 2 + g_opt.DisplayInterval = 1; + else + g_opt.DisplayInterval = 100; + end +else + g_opt.LogToConsole = 0; + g_opt.OutputFlag = 0; +end +g_opt.NonConvex = 2; +g_opt.QCPDual = 1; %% Turn on the multipliers for quadratic constraints + +if ~issparse(A) + A = sparse(A); +end +if issparse(c) + c = full(c); +end + +%% split up quadratic constraints +if ~isempty(Q) + [ieq_quad, igt_quad, ilt_quad, Q_quad, B_quad, d_quad] = ... + convert_quad_constraint(Q, B, lq, uq); + %% grab some dimensions + neq_quad = length(ieq_quad); %% number of quadratic equalities + niq_quad = length(ilt_quad) + length(igt_quad); %% number of quadratic inequalities +else + Q_quad = {}; + neq_quad = 0; niq_quad = 0; +end + +%% split up linear constraints +if isempty(Q) && ~isempty(B) + [ieq_lin, igt_lin, ilt_lin, A_lin, b_lin] = convert_lin_constraint([A; B], [l; lq], [u; uq]); +else + [ieq_lin, igt_lin, ilt_lin, A_lin, b_lin] = convert_lin_constraint(A, l, u); +end +%% grab some dimensions +neq_lin = length(ieq_lin); %% number of linear equalities +niq_lin = length(ilt_lin) + length(igt_lin); %% number of linear inequalities + +%% set up model +if ~isempty(Q_quad) + m.quadcon = cell2struct([ cellfun(@(x)(0.5*x), Q_quad, 'UniformOutput', false), ... + mat2cell(reshape(B_quad', prod(size(B_quad)) ,[]), nx*ones(neq_quad+niq_quad,1)), ... + num2cell(d_quad,2), ... + cellstr(char([double('=')*ones(neq_quad,1); double('<')*ones(niq_quad,1)]))], ... + { 'Qc' , ... + 'q' , ... + 'rhs' , ... + 'sense' }, ... + 2); +end +m.A = A_lin; +m.rhs = b_lin; +m.sense = char([ double('=')*ones(1,neq_lin) double('<')*ones(1,niq_lin) ]); +m.lb = xmin; +m.ub = xmax; +m.obj = c; +if ~isempty(vtype) + m.vtype = vtype; +end +if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I' | ... + vtype == 'S' | vtype == 'N', 1)) + mi = 0; +else + mi = 1; +end + +%% Call the solver +isemptyQ = cell2mat(cellfun(@(x)(~nnz(x)), Q_quad, 'UniformOutput', false)); +if sum(isemptyQ) == (neq_quad + niq_quad) %% No quadratic terms in quadratic constraints (linear constraints) + if ~nnz(H) + lpqcqp = 'LP'; + else + lpqcqp = 'QP'; + if ~issparse(H) + H = sparse(H); + end + m.Q = 0.5 * H; + end +else + lpqcqp = 'QCQP'; + if nnz(H) + if ~issparse(H) + H = sparse(H); + end + m.Q = 0.5 * H; + end +end +if mi + lpqcqp = ['MI' lpqcqp]; +end +if verbose + alg_names = { + 'automatic', + 'primal simplex', + 'dual simplex', + 'interior point', + 'concurrent', + 'deterministic concurrent', + 'deterministic concurrent simplex' + }; + vn = gurobiver; + fprintf('Gurobi Version %s -- %s %s solver\n', ... + vn, alg_names{g_opt.Method+2}, lpqcqp); +end + +results = gurobi(m, g_opt); + +%% Check for status of the optimization run and prepare output +switch results.status + case 'LOADED' %% 1 + eflag = -1; + case 'OPTIMAL' %% 2, optimal solution found + eflag = 1; + case 'INFEASIBLE' %% 3 + eflag = -3; + case 'INF_OR_UNBD' %% 4 + eflag = -4; + case 'UNBOUNDED' %% 5 + eflag = -5; + case 'CUTOFF' %% 6 + eflag = -6; + case 'ITERATION_LIMIT' %% 7 + eflag = -7; + case 'NODE_LIMIT' %% 8 + eflag = -8; + case 'TIME_LIMIT' %% 9 + eflag = -9; + case 'SOLUTION_LIMIT' %% 10 + eflag = -10; + case 'INTERRUPTED' %% 11 + eflag = -11; + case 'NUMERIC' %% 12 + eflag = -12; + case 'SUBOPTIMAL' %% 13 + eflag = -13; + case 'INPROGRESS' %% 14 + eflag = -14; + case 'USER_OBJ_LIMIT' %% 15 + eflag = -15; + case 'WORK_LIMIT' %% 15 + eflag = -16; + case 'MEM_LIMIT' %% 16 + eflag = -17; + otherwise + eflag = 0; +end +if nargout > 3 + output = results; +end + +%% check for empty results (in case optimization failed) +if ~isfield(results, 'x') || isempty(results.x) + x = NaN(nx, 1); + lam.lower = NaN(nx, 1); + lam.upper = NaN(nx, 1); +else + x = results.x; + lam.lower = zeros(nx, 1); + lam.upper = zeros(nx, 1); +end +if ~isfield(results, 'objval') || isempty(results.objval) + f = NaN; +else + f = results.objval; +end +if ~isfield(results, 'pi') || isempty(results.pi) + pi = NaN(length(m.rhs), 1); +else + pi = results.pi; +end +if ~isfield(results, 'rc') || isempty(results.rc) + rc = NaN(nx, 1); +else + rc = results.rc; +end + +kl = find(rc > 0); %% lower bound binding +ku = find(rc < 0); %% upper bound binding +lam.lower(kl) = rc(kl); +lam.upper(ku) = -rc(ku); + +[mu_l, mu_u] = convert_constraint_multipliers( ... + -pi(1:neq_lin), -pi(neq_lin+(1:niq_lin)), ieq_lin, igt_lin, ilt_lin); + +if ~isempty(Q_quad) + if ~isfield(results, 'qcpi') || isempty(results.qcpi) + qcpi = NaN(length(m.quadcon), 1); + % Current version of Gurobi (11.0.3) does not return multipliers for + % non-convex qcqp. See + % https://docs.gurobi.com/projects/optimizer/en/current/reference/attributes/constraintquadratic.html#qcpi + else + qcpi = results.qcpi; + end + [mu_lq, mu_uq] = convert_constraint_multipliers( ... + -qcpi(1:neq_quad), -qcpi(neq_quad+(1:niq_quad)), ... + ieq_quad, igt_quad, ilt_quad); + + lambda = struct( ... + 'mu_l' , mu_l, ... + 'mu_u' , mu_u, ... + 'mu_lq' , mu_lq, ... + 'mu_uq' , mu_uq, ... + 'lower' , lam.lower, ... + 'upper' , lam.upper ... + ); +else + lambda = struct( ... + 'mu_l' , mu_l, ... + 'mu_u' , mu_u, ... + 'lower' , lam.lower, ... + 'upper' , lam.upper ... + ); +end From a2faff6c26d1e4bfe994a74f04f1563133dea1f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilson=20Gonz=C3=A1lez=20Vanegas?= Date: Fri, 31 Jul 2026 18:00:25 -0500 Subject: [PATCH 3/4] Enable MIQCQP handling for opt_model() and mp.opt_model(). --- lib/+mp/opt_model.m | 5 ++--- lib/@opt_model/solve.m | 7 +++---- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/lib/+mp/opt_model.m b/lib/+mp/opt_model.m index adc3b30..4ada5ae 100644 --- a/lib/+mp/opt_model.m +++ b/lib/+mp/opt_model.m @@ -802,9 +802,8 @@ function clear_cached_params(obj) [x, f, eflag, output, lambda] = ... miqps_master(HH, CC, A, l, u, xmin, xmax, x0, vtype, opt); else %% MIQCQP - mixed integer quadratically constrained quadratic program - % To be implemented ... - % [x, f, eflag, output, lambda] = ... - % miqcqps_master(HH, CC, Q, B, k, ll, uu, A, l, u, xmin, xmax, x0, vtype, opt); + [x, f, eflag, output, lambda] = ... + miqcqps_master(HH, CC, Q, B, ll, uu, A, l, u, xmin, xmax, x0, vtype, opt); end else %% LP, QP - linear/quadratic program %% run solver diff --git a/lib/@opt_model/solve.m b/lib/@opt_model/solve.m index eeb75ad..46c682d 100644 --- a/lib/@opt_model/solve.m +++ b/lib/@opt_model/solve.m @@ -224,10 +224,9 @@ if isempty(Q) %% MILP, MIQP - mixed integer linear/quadratic program [x, f, eflag, output, lambda] = ... miqps_master(HH, CC, A, l, u, xmin, xmax, x0, vtype, opt); - else %% MIQCQP - mixed integer quadratically constrained quadratic program - % To be implemented ... - % [x, f, eflag, output, lambda] = ... - % miqcqps_master(HH, CC, Q, B, k, ll, uu, A, l, u, xmin, xmax, x0, vtype, opt); + else %% MIQCQP - mixed integer quadratically constrained quadratic program + [x, f, eflag, output, lambda] = ... + miqcqps_master(HH, CC, Q, B, ll, uu, A, l, u, xmin, xmax, x0, vtype, opt); end else %% LP, QP - linear/quadratic program %% run solver From 935148436ec1ca2b4f692c42455387b6cd71e167 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilson=20Gonz=C3=A1lez=20Vanegas?= Date: Fri, 31 Jul 2026 18:02:47 -0500 Subject: [PATCH 4/4] Add test scripts for MIQP/MIQCQP small problems. --- lib/t/t_miqcqps_master.m | 198 ++++++++++++++++++++++++++++++++ lib/t/t_mm_solve_miqcqps.m | 221 ++++++++++++++++++++++++++++++++++++ lib/t/t_om_solve_miqcqps.m | 225 +++++++++++++++++++++++++++++++++++++ 3 files changed, 644 insertions(+) create mode 100644 lib/t/t_miqcqps_master.m create mode 100644 lib/t/t_mm_solve_miqcqps.m create mode 100644 lib/t/t_om_solve_miqcqps.m diff --git a/lib/t/t_miqcqps_master.m b/lib/t/t_miqcqps_master.m new file mode 100644 index 0000000..2e20eb7 --- /dev/null +++ b/lib/t/t_miqcqps_master.m @@ -0,0 +1,198 @@ +function t_miqcqps_master(quiet) +% t_miqcqps_master - Tests of MIQP/MIQCQP solvers via miqcqps_master. + +% MP-Opt-Model +% Copyright (c) 2019-2026, Power Systems Engineering Research Center (PSERC) +% by Wilson Gonzalez Vanegas, Universidad Nacional de Colombia Sede Manizales +% and Ray Zimmerman, PSERC Cornell +% +% This file is part of MP-Opt-Model. +% Covered by the 3-clause BSD License (see LICENSE file for details). +% See https://github.com/MATPOWER/mp-opt-model for more info. + +if nargin < 1 + quiet = 0; +end + +% 1 2 3 4 +algs = {'DEFAULT', 'CPLEX', 'GUROBI', 'MOSEK'}; +names = {'DEFAULT', 'CPLEX', 'Gurobi', 'MOSEK'}; +check = { [], 'cplex', 'gurobi', 'mosek'}; + +% 1 2 3 4 +does_miqp = [1 1 1 1]; +does_miqcqp = [1 1 0 0]; +does_nonconv = [1 1 0 0]; + +nmiqp_convex = 3; +nmiqp_nonconvex = 6; +nmiqcqp_convex = 6; +nmiqcqp_nonconvex = 9; + +n = nmiqp_convex+nmiqp_nonconvex+nmiqcqp_convex+nmiqcqp_nonconvex; + +t_begin(n * length(algs), quiet); + +for k = 1:length(algs) + if ~isempty(check{k}) && ~have_feature(check{k}) + t_skip(n, sprintf('%s not installed', names{k})); + else + opt = struct('verbose', 0, 'alg', algs{k}); + + mpopt = struct( ... + 'verbose', 0, ... + 'opf', struct( ... + 'violation', 1e-6 ), ... + 'gurobi', struct( ... + 'method', -1, ... + 'timelimit', Inf, ... + 'threads', 0, ... + 'opts', [], ... + 'opt_fname', [], ... + 'opt', 0) ... + ); + if have_feature('gurobi') + opt.grb_opt = gurobi_options([], mpopt); + opt.grb_opt.BarQCPConvTol = 1e-8; + opt.gurobi.MIPGap = 1e-9; + opt.gurobi.MIPGapAbs = 1e-9; + opt.gurobi.FeasibilityTol = 1e-9; + opt.gurobi.OptimalityTol = 1e-9; + opt.gurobi.IntFeasTol = 1e-9; + end + opt_r = opt; + opt_r.relax_integer = 1; + opt_f = opt; + opt_f.fix_integer = 1; + + if does_miqp(k) + t = sprintf('%s - convex 2-d MIQP : ', names{k}); + %% 1) From OPTI Toolbox: https://jonathancurrie.github.io/OPTI/examples/problem-types/miqp/ + H = [1 -1; -1 2]; + c = [-2 -6]'; + Q = {}; B = []; lq = []; uq = []; + A = [1 1; -1 2; 2 1]; + u = [2 2 3]'; + l = []; + xmin = zeros(2,1); + xmax = Inf(2,1); + x0 = []; + vtype = 'IC'; + [x, f, s, out, lam] = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt); + t_is(s, 1, 12, [t 'success']); + t_is(x, [1, 1]', 7, [t 'x']); + t_is(f, -7.500, 4, [t 'f']); + + %% 2) From Billionnet et al. 2012 (https://hal.science/hal-01125718v1/file/Extending_the_QCR_method_to_general_mixed-integer_.pdf) + t = sprintf('%s - nonconvex 4-d MIQP : ', names{k}); + if does_nonconv(k) + H = [ -14 6 -30 -8; + 6 -28 -14 -26; + -30 -14 16 14; + -8 -26 14 24 ]; + c = [15; 10; -7; -4]; + Q = {}; B = []; lq = []; uq = []; + A = [5 1 8 4]; + u = 95; + l = []; + xmin = zeros(4,1); + xmax = 10*ones(4,1); + x0 = []; + vtype = 'IICC'; + [x, f, s, out, lam] = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt); + t_is(s, 1, 12, [t 'success']); + t_is(x, [8, 10, 2.0268, 7.1964]', 2, [t 'x']); + t_is(f, -3434.2700, 2, [t 'f']); + + %% 3) Same previous fixing integer variables to their optimal values + t = sprintf('%s - nonconvex 4-d MIQP (integer fixed) : ', names{k}); + x0 = [8, 10, 0, 0]'; + [x, f, s, out, lam] = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt_f); + t_is(s, 1, 12, [t 'success']); + t_is(x, [8, 10, 2.0268, 7.1964]', 2, [t 'x']); + t_is(f, -3434.2700, 2, [t 'f']); + else + t_skip(nmiqp_nonconvex, sprintf('%s does not handle nonconvex MIQP problems', names{k})); + end + else + t_skip(nmiqp_convex+nmiqp_nonconvex, sprintf('%s does not handle MIQP problems', names{k})); + end + + %% 4) From OPTI Toolbox: https://jonathancurrie.github.io/OPTI/examples/problem-types/miqcqp/ + if does_miqcqp(k) + t = sprintf('%s - convex 2-d MIQCQP : ', names{k}); + H = [1 0; 0 1]; + c = [-2 -2]'; + Q = {2*speye(2)}; + B = [0 -2]; + lq = []; + uq = 1; + A = [-1 1; 1 3]; + u = [2 5]'; + l = []; + xmin = zeros(2,1); + xmax = Inf(2,1); + x0 = []; + vtype = 'IC'; + [x, f, s, out, lam] = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt); + t_is(s, 1, 12, [t 'success']); + t_is(x, [1.0, 4/3]', 7, [t 'x']); + t_is(f, -3.2778, 4, [t 'f']); + + %% 5) Same previous relaxing integer variable + t = sprintf('%s - convex 2-d MIQCQP (integer relaxed) : ', names{k}); + [x, f, s, out, lam] = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt_r); + t_is(s, 1, 12, [t 'success']); + t_is(x, [7/5, 6/5]', 7, [t 'x']); + t_is(f, -3.5000, 4, [t 'f']); + + %% 6) From You and Dai 2020 (https://ieeexplore.ieee.org/abstract/document/9304368) + %% [31/07/2026] - WGV: the sign of -24.3 for variable x3 in second row of A was corrected + %% after reviewing the main paper of 1960 mentioned in You and Dai. The gold + %% standard results are taken also from the paper of 1960 (Land and Doig) + t = sprintf('%s - nonconvex 7-d MIQCQP : ', names{k}); + if does_nonconv(k) + H = sparse([1;2],[2;1],[1/2;1/2],7,7,2); + c = [-77.9 -76.8 -89.6 -97.1 -31.3 0 0]'; + Q = {sparse([1;2;3;3],[3;3;1;2],[-1/2;1/2;-1/2;1/2],7,7,4)}; + B = zeros(1,7); + lq = []; + uq = 0; + A = [10.9 3.6 -40.8 43.9 7.1 1 0 + -86.8 32.7 24.3 13.8 -12.6 0 1 + 60.9 68.9 69.0 -56.9 22.5 0 0]; + u = [82.3 77.3 86.5]'; + l = u; + xmin = zeros(7,1); + xmax = Inf(7,1); + x0 = []; + vtype = 'IIICCCC'; + [x, f, s, out, lam] = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt); + t_is(s, 1, 12, [t 'success']); + t_is(x([1:6]), [1 0 4 5.0702 1.6930 0]', 4, [t 'x']); + t_is(f, -981.6023, 4, [t 'f']); + + %% 7) Same previous relaxing integer variable + t = sprintf('%s - nonconvex 7-d MIQCQP (integer relaxed): ', names{k}); + [x, f, s, out, lam] = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt_r); + t_is(s, 1, 12, [t 'success']); + t_is(x, [1.4960 0 5.0210 6.1697 0 0 0]', 4, [t 'x']); + t_is(f, -1165.50, 2, [t 'f']); + + %% 8) Same problem fixing integer variables to their optimal values + t = sprintf('%s - nonconvex 7-d MIQCQP (integer fixed) : ', names{k}); + x0 = [1 0 4 zeros(1,4)]'; + [x, f, s, out, lam] = miqcqps_master(H, c, Q, B, lq, uq, A, l, u, xmin, xmax, x0, vtype, opt_f); + t_is(s, 1, 12, [t 'success']); + t_is(x([1:6]), [1 0 4 5.0702 1.6930 0]', 4, [t 'x']); + t_is(f, -981.6023, 4, [t 'f']); + else + t_skip(nmiqcqp_nonconvex, sprintf('%s does not handle nonconvex MIQCQP problems', names{k})); + end + else + t_skip(nmiqcqp_convex+nmiqcqp_nonconvex, sprintf('%s does not handle MIQCQP problems', names{k})); + end + end +end + +t_end; diff --git a/lib/t/t_mm_solve_miqcqps.m b/lib/t/t_mm_solve_miqcqps.m new file mode 100644 index 0000000..0c4a6c7 --- /dev/null +++ b/lib/t/t_mm_solve_miqcqps.m @@ -0,0 +1,221 @@ +function t_mm_solve_miqcqps(quiet) +% t_mm_solve_miqcqps - Tests of MIQP/MIQCQP solvers via mp.opt_model.solve. + +% MP-Opt-Model +% Copyright (c) 2019-2026, Power Systems Engineering Research Center (PSERC) +% by Wilson Gonzalez Vanegas, Universidad Nacional de Colombia Sede Manizales +% and Ray Zimmerman, PSERC Cornell +% +% This file is part of MP-Opt-Model. +% Covered by the 3-clause BSD License (see LICENSE file for details). +% See https://github.com/MATPOWER/mp-opt-model for more info. + +if nargin < 1 + quiet = 0; +end + +% 1 2 3 4 +algs = {'DEFAULT', 'CPLEX', 'GUROBI', 'MOSEK'}; +names = {'DEFAULT', 'CPLEX', 'Gurobi', 'MOSEK'}; +check = { [], 'cplex', 'gurobi', 'mosek'}; + +% 1 2 3 4 +does_miqp = [1 1 1 1]; +does_miqcqp = [1 1 0 0]; +does_nonconv = [1 1 0 0]; + +nmiqp_convex = 3; +nmiqp_nonconvex = 6; +nmiqcqp_convex = 6; +nmiqcqp_nonconvex = 9; + +n = nmiqp_convex+nmiqp_nonconvex+nmiqcqp_convex+nmiqcqp_nonconvex; + +t_begin(n * length(algs), quiet); + +for k = 1:length(algs) + if ~isempty(check{k}) && ~have_feature(check{k}) + t_skip(n, sprintf('%s not installed', names{k})); + else + opt = struct('verbose', 0, 'alg', algs{k}); + + mpopt = struct( ... + 'verbose', 0, ... + 'opf', struct( ... + 'violation', 1e-6 ), ... + 'gurobi', struct( ... + 'method', -1, ... + 'timelimit', Inf, ... + 'threads', 0, ... + 'opts', [], ... + 'opt_fname', [], ... + 'opt', 0) ... + ); + if have_feature('gurobi') + opt.grb_opt = gurobi_options([], mpopt); + opt.grb_opt.BarQCPConvTol = 1e-8; + opt.gurobi.MIPGap = 1e-9; + opt.gurobi.MIPGapAbs = 1e-9; + opt.gurobi.FeasibilityTol = 1e-9; + opt.gurobi.OptimalityTol = 1e-9; + opt.gurobi.IntFeasTol = 1e-9; + end + opt_r = opt; + opt_r.relax_integer = 1; + opt_f = opt; + opt_f.fix_integer = 1; + + if does_miqp(k) + t = sprintf('%s - convex 2-d MIQP : ', names{k}); + %% 1) From OPTI Toolbox: https://jonathancurrie.github.io/OPTI/examples/problem-types/miqp/ + H = [1 -1; -1 2]; + c = [-2 -6]'; + Q = {spalloc(2,2,0)}; B = []; lq = []; uq = []; + A = [1 1; -1 2; 2 1]; + u = [2 2 3]'; + l = []; + xmin = zeros(2,1); + xmax = Inf(2,1); + x0 = []; + vtype = 'IC'; + mm = mp.opt_model(); + mm.var.add('x', 2, x0, xmin, xmax, vtype); + mm.qdc.add(mm.var, 'cost', H, c); + mm.qcn.add(mm.var, 'xQx', Q, B, lq, uq); + mm.lin.add(mm.var, 'Ax', A, l, u); + [x, f, s, out, lam] = mm.solve(opt); + t_is(s, 1, 12, [t 'success']); + t_is(x, [1, 1]', 7, [t 'x']); + t_is(f, -7.500, 4, [t 'f']); + + %% 2) From Billionnet et al. 2012 (https://hal.science/hal-01125718v1/file/Extending_the_QCR_method_to_general_mixed-integer_.pdf) + t = sprintf('%s - nonconvex 4-d MIQP : ', names{k}); + if does_nonconv(k) + H = [ -14 6 -30 -8; + 6 -28 -14 -26; + -30 -14 16 14; + -8 -26 14 24 ]; + c = [15; 10; -7; -4]; + Q = {spalloc(4,4,0)}; B = []; lq = []; uq = []; + A = [5 1 8 4]; + u = 95; + l = []; + xmin = zeros(4,1); + xmax = 10*ones(4,1); + x0 = []; + vtype = 'IICC'; + mm = mp.opt_model(); + mm.var.add('x', 4, x0, xmin, xmax, vtype); + mm.qdc.add(mm.var, 'cost', H, c); + mm.qcn.add(mm.var, 'xQx', Q, B, lq, uq); + mm.lin.add(mm.var, 'Ax', A, l, u); + [x, f, s, out, lam] = mm.solve(opt); + t_is(s, 1, 12, [t 'success']); + t_is(x, [8, 10, 2.0268, 7.1964]', 2, [t 'x']); + t_is(f, -3434.2700, 2, [t 'f']); + + %% 3) Same previous fixing integer variables to their optimal values + t = sprintf('%s - nonconvex 4-d MIQP (integer fixed) : ', names{k}); + x0 = [8, 10, 0, 0]'; + mm.var.set_params('x','v0',x0); + [x, f, s, out, lam] = mm.solve(opt_f); + t_is(s, 1, 12, [t 'success']); + t_is(x, [8, 10, 2.0268, 7.1964]', 2, [t 'x']); + t_is(f, -3434.2700, 2, [t 'f']); + else + t_skip(nmiqp_nonconvex, sprintf('%s does not handle nonconvex MIQP problems', names{k})); + end + else + t_skip(nmiqp_convex+nmiqp_nonconvex, sprintf('%s does not handle MIQP problems', names{k})); + end + + %% 4) From OPTI Toolbox: https://jonathancurrie.github.io/OPTI/examples/problem-types/miqcqp/ + if does_miqcqp(k) + t = sprintf('%s - convex 2-d MIQCQP : ', names{k}); + H = [1 0; 0 1]; + c = [-2 -2]'; + Q = {2*speye(2)}; + B = [0 -2]; + lq = []; + uq = 1; + A = [-1 1; 1 3]; + u = [2 5]'; + l = []; + xmin = zeros(2,1); + xmax = Inf(2,1); + x0 = []; + vtype = 'IC'; + mm = mp.opt_model(); + mm.var.add('x', 2, x0, xmin, xmax, vtype); + mm.qdc.add(mm.var, 'cost', H, c); + mm.qcn.add(mm.var, 'xQx', Q, B, lq, uq); + mm.lin.add(mm.var, 'Ax', A, l, u); + [x, f, s, out, lam] = mm.solve(opt); + t_is(s, 1, 12, [t 'success']); + t_is(x, [1.0, 4/3]', 7, [t 'x']); + t_is(f, -3.2778, 4, [t 'f']); + + %% 5) Same previous relaxing integer variable + t = sprintf('%s - convex 2-d MIQCQP (integer relaxed) : ', names{k}); + [x, f, s, out, lam] = mm.solve(opt_r); + t_is(s, 1, 12, [t 'success']); + t_is(x, [7/5, 6/5]', 7, [t 'x']); + t_is(f, -3.5000, 4, [t 'f']); + + %% 6) From You and Dai 2020 (https://ieeexplore.ieee.org/abstract/document/9304368) + %% [31/07/2026] - WGV: the sign of -24.3 for variable x3 in second row of A was corrected. + %% After reviewing, the main paper of 1960 mentioned in You and Dai. The gold + %% standard results are taken also from the paper of 1960 (Land and Doig) + t = sprintf('%s - nonconvex 7-d MIQCQP : ', names{k}); + if does_nonconv(k) + H = sparse([1;2],[2;1],[1/2;1/2],7,7,2); + c = [-77.9 -76.8 -89.6 -97.1 -31.3 0 0]'; + Q = {sparse([1;2;3;3],[3;3;1;2],[-1/2;1/2;-1/2;1/2],7,7,4)}; + B = zeros(1,7); + lq = []; + uq = 0; + A = [10.9 3.6 -40.8 43.9 7.1 1 0 + -86.8 32.7 24.3 13.8 -12.6 0 1 + 60.9 68.9 69.0 -56.9 22.5 0 0]; + u = [82.3 77.3 86.5]'; + l = u; + xmin = zeros(7,1); + xmax = Inf(7,1); + x0 = []; + vtype = 'IIICCCC'; + mm = mp.opt_model(); + mm.var.add('x', 7, x0, xmin, xmax, vtype); + mm.qdc.add(mm.var, 'cost', H, c); + mm.qcn.add(mm.var, 'xQx', Q, B, lq, uq); + mm.lin.add(mm.var, 'Ax', A, l, u); + [x, f, s, out, lam] = mm.solve(opt); + t_is(s, 1, 12, [t 'success']); + t_is(x([1:6]), [1 0 4 5.0702 1.6930 0]', 4, [t 'x']); + t_is(f, -981.6023, 4, [t 'f']); + + %% 7) Same previous relaxing integer variable + t = sprintf('%s - nonconvex 7-d MIQCQP (integer relaxed): ', names{k}); + mm.var.set_params('x','v0',ones(7,1)); + [x, f, s, out, lam] = mm.solve(opt_r); + t_is(s, 1, 12, [t 'success']); + t_is(x, [1.4960 0 5.0210 6.1697 0 0 0]', 4, [t 'x']); + t_is(f, -1165.50, 2, [t 'f']); + + %% 8) Same problem fixing integer variables to their optimal values + t = sprintf('%s - nonconvex 7-d MIQCQP (integer fixed) : ', names{k}); + x0 = [1 0 4 zeros(1,4)]'; + mm.var.set_params('x','v0',x0); + [x, f, s, out, lam] = mm.solve(opt_f); + t_is(s, 1, 12, [t 'success']); + t_is(x([1:6]), [1 0 4 5.0702 1.6930 0]', 4, [t 'x']); + t_is(f, -981.6023, 4, [t 'f']); + else + t_skip(nmiqcqp_nonconvex, sprintf('%s does not handle nonconvex MIQCQP problems', names{k})); + end + else + t_skip(nmiqcqp_convex+nmiqcqp_nonconvex, sprintf('%s does not handle MIQCQP problems', names{k})); + end + end +end + +t_end; diff --git a/lib/t/t_om_solve_miqcqps.m b/lib/t/t_om_solve_miqcqps.m new file mode 100644 index 0000000..c45289f --- /dev/null +++ b/lib/t/t_om_solve_miqcqps.m @@ -0,0 +1,225 @@ +function t_om_solve_miqcqps(quiet) +% t_om_solve_miqcqps - Tests of MIQP/MIQCQP solvers via opt_model.solve. + +% MP-Opt-Model +% Copyright (c) 2019-2026, Power Systems Engineering Research Center (PSERC) +% by Wilson Gonzalez Vanegas, Universidad Nacional de Colombia Sede Manizales +% and Ray Zimmerman, PSERC Cornell +% +% This file is part of MP-Opt-Model. +% Covered by the 3-clause BSD License (see LICENSE file for details). +% See https://github.com/MATPOWER/mp-opt-model for more info. + +if nargin < 1 + quiet = 0; +end + +% 1 2 3 4 +algs = {'DEFAULT', 'CPLEX', 'GUROBI', 'MOSEK'}; +names = {'DEFAULT', 'CPLEX', 'Gurobi', 'MOSEK'}; +check = { [], 'cplex', 'gurobi', 'mosek'}; + +% 1 2 3 4 +does_miqp = [1 1 1 1]; +does_miqcqp = [1 1 0 0]; +does_nonconv = [1 1 0 0]; + +nmiqp_convex = 3; +nmiqp_nonconvex = 6; +nmiqcqp_convex = 6; +nmiqcqp_nonconvex = 9; + +n = nmiqp_convex+nmiqp_nonconvex+nmiqcqp_convex+nmiqcqp_nonconvex; + +t_begin(n * length(algs), quiet); + +for k = 1:length(algs) + if ~isempty(check{k}) && ~have_feature(check{k}) + t_skip(n, sprintf('%s not installed', names{k})); + else + opt = struct('verbose', 0, 'alg', algs{k}); + + mpopt = struct( ... + 'verbose', 0, ... + 'opf', struct( ... + 'violation', 1e-6 ), ... + 'gurobi', struct( ... + 'method', -1, ... + 'timelimit', Inf, ... + 'threads', 0, ... + 'opts', [], ... + 'opt_fname', [], ... + 'opt', 0) ... + ); + if have_feature('gurobi') + opt.grb_opt = gurobi_options([], mpopt); + opt.grb_opt.BarQCPConvTol = 1e-8; + opt.gurobi.MIPGap = 1e-9; + opt.gurobi.MIPGapAbs = 1e-9; + opt.gurobi.FeasibilityTol = 1e-9; + opt.gurobi.OptimalityTol = 1e-9; + opt.gurobi.IntFeasTol = 1e-9; + end + opt_r = opt; + opt_r.relax_integer = 1; + opt_f = opt; + opt_f.fix_integer = 1; + + if does_miqp(k) + t = sprintf('%s - convex 2-d MIQP : ', names{k}); + %% 1) From OPTI Toolbox: https://jonathancurrie.github.io/OPTI/examples/problem-types/miqp/ + H = [1 -1; -1 2]; + c = [-2 -6]'; + Q = {spalloc(2,2,0)}; B = []; lq = []; uq = []; + A = [1 1; -1 2; 2 1]; + u = [2 2 3]'; + l = []; + xmin = zeros(2,1); + xmax = Inf(2,1); + x0 = []; + vtype = 'IC'; + om = opt_model(); + om.init_set_types(); + om.var.add('x', 2, x0, xmin, xmax, vtype); + om.qdc.add(om.var, 'cost', H, c); + om.qcn.add(om.var, 'xQx', Q, B, lq, uq); + om.lin.add(om.var, 'Ax', A, l, u); + [x, f, s, out, lam] = om.solve(opt); + t_is(s, 1, 12, [t 'success']); + t_is(x, [1, 1]', 7, [t 'x']); + t_is(f, -7.500, 4, [t 'f']); + + %% 2) From Billionnet et al. 2012 (https://hal.science/hal-01125718v1/file/Extending_the_QCR_method_to_general_mixed-integer_.pdf) + t = sprintf('%s - nonconvex 4-d MIQP : ', names{k}); + if does_nonconv(k) + H = [ -14 6 -30 -8; + 6 -28 -14 -26; + -30 -14 16 14; + -8 -26 14 24 ]; + c = [15; 10; -7; -4]; + Q = {spalloc(4,4,0)}; B = []; lq = []; uq = []; + A = [5 1 8 4]; + u = 95; + l = []; + xmin = zeros(4,1); + xmax = 10*ones(4,1); + x0 = []; + vtype = 'IICC'; + om = opt_model(); + om.init_set_types(); + om.var.add('x', 4, x0, xmin, xmax, vtype); + om.qdc.add(om.var, 'cost', H, c); + om.qcn.add(om.var, 'xQx', Q, B, lq, uq); + om.lin.add(om.var, 'Ax', A, l, u); + [x, f, s, out, lam] = om.solve(opt); + t_is(s, 1, 12, [t 'success']); + t_is(x, [8, 10, 2.0268, 7.1964]', 2, [t 'x']); + t_is(f, -3434.2700, 2, [t 'f']); + + %% 3) Same previous fixing integer variables to their optimal values + t = sprintf('%s - nonconvex 4-d MIQP (integer fixed) : ', names{k}); + x0 = [8, 10, 0, 0]'; + om.var.set_params('x','v0',x0); + [x, f, s, out, lam] = om.solve(opt_f); + t_is(s, 1, 12, [t 'success']); + t_is(x, [8, 10, 2.0268, 7.1964]', 2, [t 'x']); + t_is(f, -3434.2700, 2, [t 'f']); + else + t_skip(nmiqp_nonconvex, sprintf('%s does not handle nonconvex MIQP problems', names{k})); + end + else + t_skip(nmiqp_convex+nmiqp_nonconvex, sprintf('%s does not handle MIQP problems', names{k})); + end + + %% 4) From OPTI Toolbox: https://jonathancurrie.github.io/OPTI/examples/problem-types/miqcqp/ + if does_miqcqp(k) + t = sprintf('%s - convex 2-d MIQCQP : ', names{k}); + H = [1 0; 0 1]; + c = [-2 -2]'; + Q = {2*speye(2)}; + B = [0 -2]; + lq = []; + uq = 1; + A = [-1 1; 1 3]; + u = [2 5]'; + l = []; + xmin = zeros(2,1); + xmax = Inf(2,1); + x0 = []; + vtype = 'IC'; + om = opt_model(); + om.init_set_types(); + om.var.add('x', 2, x0, xmin, xmax, vtype); + om.qdc.add(om.var, 'cost', H, c); + om.qcn.add(om.var, 'xQx', Q, B, lq, uq); + om.lin.add(om.var, 'Ax', A, l, u); + [x, f, s, out, lam] = om.solve(opt); + t_is(s, 1, 12, [t 'success']); + t_is(x, [1.0, 4/3]', 7, [t 'x']); + t_is(f, -3.2778, 4, [t 'f']); + + %% 5) Same previous relaxing integer variable + t = sprintf('%s - convex 2-d MIQCQP (integer relaxed) : ', names{k}); + [x, f, s, out, lam] = om.solve(opt_r); + t_is(s, 1, 12, [t 'success']); + t_is(x, [7/5, 6/5]', 7, [t 'x']); + t_is(f, -3.5000, 4, [t 'f']); + + %% 6) From You and Dai 2020 (https://ieeexplore.ieee.org/abstract/document/9304368) + %% [31/07/2026] - WGV: the sign of -24.3 for variable x3 in second row of A was corrected. + %% After reviewing, the main paper of 1960 mentioned in You and Dai. The gold + %% standard results are taken also from the paper of 1960 (Land and Doig) + t = sprintf('%s - nonconvex 7-d MIQCQP : ', names{k}); + if does_nonconv(k) + H = sparse([1;2],[2;1],[1/2;1/2],7,7,2); + c = [-77.9 -76.8 -89.6 -97.1 -31.3 0 0]'; + Q = {sparse([1;2;3;3],[3;3;1;2],[-1/2;1/2;-1/2;1/2],7,7,4)}; + B = zeros(1,7); + lq = []; + uq = 0; + A = [10.9 3.6 -40.8 43.9 7.1 1 0 + -86.8 32.7 24.3 13.8 -12.6 0 1 + 60.9 68.9 69.0 -56.9 22.5 0 0]; + u = [82.3 77.3 86.5]'; + l = u; + xmin = zeros(7,1); + xmax = Inf(7,1); + x0 = []; + vtype = 'IIICCCC'; + om = opt_model(); + om.init_set_types(); + om.var.add('x', 7, x0, xmin, xmax, vtype); + om.qdc.add(om.var, 'cost', H, c); + om.qcn.add(om.var, 'xQx', Q, B, lq, uq); + om.lin.add(om.var, 'Ax', A, l, u); + [x, f, s, out, lam] = om.solve(opt); + t_is(s, 1, 12, [t 'success']); + t_is(x([1:6]), [1 0 4 5.0702 1.6930 0]', 4, [t 'x']); + t_is(f, -981.6023, 4, [t 'f']); + + %% 7) Same previous relaxing integer variable + t = sprintf('%s - nonconvex 7-d MIQCQP (integer relaxed): ', names{k}); + om.var.set_params('x','v0',ones(7,1)); + [x, f, s, out, lam] = om.solve(opt_r); + t_is(s, 1, 12, [t 'success']); + t_is(x, [1.4960 0 5.0210 6.1697 0 0 0]', 4, [t 'x']); + t_is(f, -1165.50, 2, [t 'f']); + + %% 8) Same problem fixing integer variables to their optimal values + t = sprintf('%s - nonconvex 7-d MIQCQP (integer fixed) : ', names{k}); + x0 = [1 0 4 zeros(1,4)]'; + om.var.set_params('x','v0',x0); + [x, f, s, out, lam] = om.solve(opt_f); + t_is(s, 1, 12, [t 'success']); + t_is(x([1:6]), [1 0 4 5.0702 1.6930 0]', 4, [t 'x']); + t_is(f, -981.6023, 4, [t 'f']); + else + t_skip(nmiqcqp_nonconvex, sprintf('%s does not handle nonconvex MIQCQP problems', names{k})); + end + else + t_skip(nmiqcqp_convex+nmiqcqp_nonconvex, sprintf('%s does not handle MIQCQP problems', names{k})); + end + end +end + +t_end;