Documentation

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

Suppose that you can put your problem in 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, and the Optimization Decision Table states that `fmincon` is the appropriate solver.

The example assumes that the quadratic matrices are symmetric. This is 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 the Hi are N-by-N matrices, f and the ki are N-by-1 vectors, and c and the 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.)

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

For 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. Examine 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` used just 10 function evaluations to solve the problem.

Compare this 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...'```

This time `fmincon` used 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 Ji, pi, and qi variables. The `lambda.eqnonlin(i)` structure has the Lagrange multipliers for equality constraints.