Documentation

Quadratic Constrained Problem

Suppose that your problem has the form

subject to

where 1 ≤ i ≤ m. Assume that at least one H_i is nonzero; otherwise, you can use quadprog or linprog to solve this problem. With nonzero H_i, the constraints are nonlinear, which means fmincon is the appropriate solver according to the Optimization Decision Table.

The example assumes that the quadratic matrices are symmetric without loss of generality. You can replace a nonsymmetric H (or Q) matrix with an equivalent symmetrized version (H + H^T)/2.

If x has N components, then Q and H_i are N-by-N matrices, f and k_i are N-by-1 vectors, and c and d_i are scalars.

Objective Function

Formulate the problem using fmincon syntax. Assume that x and f are column vectors. (x is a column vector when the initial vector x0 is a column vector.)

function [y,grady] = quadobj(x,Q,f,c)
y = 1/2*x'*Q*x + f'*x + c;
if nargout > 1
    grady = Q*x + f;
end

Constraint Function

For consistency and easy indexing, place every quadratic constraint matrix in one cell array. Similarly, place the linear and constant terms in cell arrays.

function [y,yeq,grady,gradyeq] = quadconstr(x,H,k,d)
jj = length(H); % jj is the number of inequality constraints
y = zeros(1,jj);
for i = 1:jj
    y(i) = 1/2*x'*H{i}*x + k{i}'*x + d{i};
end
yeq = [];
    
if nargout > 2
    grady = zeros(length(x),jj);
    for i = 1:jj
        grady(:,i) = H{i}*x + k{i};
    end
end
gradyeq = [];

Numeric Example

Suppose that you have the following problem.

Q = [3,2,1;
     2,4,0;
     1,0,5];
f = [-24;-48;-130];
c = -2;

rng default % For reproducibility
% Two sets of random quadratic constraints:
H{1} = gallery('randcorr',3); % Random positive definite matrix
H{2} = gallery('randcorr',3);
k{1} = randn(3,1);
k{2} = randn(3,1);
d{1} = randn;
d{2} = randn;

Hessian

Create a Hessian function. The Hessian of the Lagrangian is given by the equation

fmincon calculates an approximate set of Lagrange multipliers λ_i, and packages them in a structure. To include the Hessian, use the following function.

function hess = quadhess(x,lambda,Q,H)
hess = Q;
jj = length(H); % jj is the number of inequality constraints
for i = 1:jj
    hess = hess + lambda.ineqnonlin(i)*H{i};
end

Solution

Use the fmincon interior-point algorithm to solve the problem most efficiently. This algorithm accepts a Hessian function that you supply. Set these options.

options = optimoptions(@fmincon,'Algorithm','interior-point',...
    'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,...
    'HessianFcn',@(x,lambda)quadhess(x,lambda,Q,H));

Call fmincon to solve the problem.

fun = @(x)quadobj(x,Q,f,c);
nonlconstr = @(x)quadconstr(x,H,k,d);
x0 = [0;0;0]; % Column vector
[x,fval,eflag,output,lambda] = fmincon(fun,x0,...
    [],[],[],[],[],[],nonlconstr,options);

Examine the Lagrange multipliers.

lambda.ineqnonlin

ans =

   12.8412
   39.2337

Both nonlinear inequality multipliers are nonzero, so both quadratic constraints are active at the solution.

Efficiency When Providing a Hessian

The interior-point algorithm with gradients and a Hessian is efficient. View the number of function evaluations.

output

output = 

         iterations: 9
          funcCount: 10
    constrviolation: 0
           stepsize: 5.3547e-04
          algorithm: 'interior-point'
      firstorderopt: 1.5851e-05
       cgiterations: 0
            message: 'Local minimum found that satisfies the constraints.

Optimization compl...'

fmincon takes only 10 function evaluations to solve the problem.

Compare this result to the solution without the Hessian.

options.HessianFcn = [];
[x2,fval2,eflag2,output2,lambda2] = fmincon(fun,[0;0;0],...
    [],[],[],[],[],[],nonlconstr,options);
output2

output2 = 

         iterations: 17
          funcCount: 22
    constrviolation: 0
           stepsize: 2.8475e-04
          algorithm: 'interior-point'
      firstorderopt: 1.7680e-05
       cgiterations: 0
            message: 'Local minimum found that satisfies the constraints.

Optimization compl...'

In this case, fmincon takes about twice as many iterations and function evaluations. The solutions are the same to within tolerances.

Extension to Quadratic Equality Constraints

If you also have quadratic equality constraints, you can use essentially the same technique. The problem is the same, with the additional constraints

Reformulate your constraints to use the J_i, p_i, and q_i variables. The lambda.eqnonlin(i) structure has the Lagrange multipliers for equality constraints.