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 Hi is nonzero; otherwise, you can use `quadprog` or `linprog` to solve this problem. With nonzero Hi, 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 + HT)/2.

If x has N components, then Q and Hi are N-by-N matrices, f and ki are N-by-1 vectors, and c and di 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

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

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

`$\frac{1}{2}{x}^{T}{J}_{i}x+{p}_{i}^{T}x+{q}_{i}=0.$`
Reformulate your constraints to use the Ji, pi, and qi variables. The `lambda.eqnonlin(i)` structure has the Lagrange multipliers for equality constraints.