Main Content

This example shows how to solve an optimization problem that has a linear or quadratic objective and quadratic inequality constraints. The example generates and uses the gradient and Hessian of the objective and constraint functions.

Suppose that your problem has the form

$$\underset{x}{\mathrm{min}}\left(\frac{1}{2}{x}^{T}Qx+{f}^{T}x+c\right)$$

subject to

$$\frac{1}{2}{x}^{T}{H}_{i}x+{k}_{i}^{T}x+{d}_{i}\le 0,$$

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 `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

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

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 = [];

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;

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

$${\nabla}_{xx}^{2}L(x,\lambda )={\nabla}^{2}f(x)+{\displaystyle \sum {\lambda}_{i}{\nabla}^{2}{c}_{i}(x)}+{\displaystyle \sum {\lambda}_{i}{\nabla}^{2}ce{q}_{i}(x)}.$$

`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

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.

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.

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

$$\frac{1}{2}{x}^{T}{J}_{i}x+{p}_{i}^{T}x+{q}_{i}=0.$$

Reformulate your constraints to use the *J _{i}*,

`lambda.eqnonlin(i)`

structure
has the Lagrange multipliers for equality constraints.