## Linear or Quadratic Objective with Quadratic Constraints

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.

### Quadratic Constrained Problem

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 *H*, the constraints are nonlinear, which means

_{i}`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*are

_{i}*N*-by-1 vectors, and

*c*and

*d*are scalars.

_{i}### 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

$${\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 puts 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

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

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

*p*, and

_{i}*q*variables. The

_{i}`lambda.eqnonlin(i)`

structure
has the Lagrange multipliers for equality constraints.