Accelerating the pace of engineering and science

# Optimization Toolbox

## Large-Scale Bound Constrained Quadratic Programming

This example shows how to determine the shape of a circus tent by solving a large-scale quadratic optimization problem. The shape of a circus tent is determined by a constrained optimization problem. We will solve this problem with the function `quadprog` in Optimization Toolbox™.

The Problem

Imagine building a circus tent to cover a square lot. The tent has five poles that will be covered with an elastic material. From this structure, we want to find the natural shape of the tent. This natural shape corresponds to the minimum of a certain energy function computed from the surface position and squared norm of its gradient.

```% Draw the tent poles. largeL = zeros(36); mask = [6 7 30 31]; largeL(mask,mask) = .3*ones(4); largeL(18:19,18:19) = .5*ones(2); xx = [1:5,5:6,6:15,15:16,16:25,25:26,26:30]; [XX,YY] = meshgrid(xx) ; axis([1 30 1 30 0 .5],'off'); surface(XX,YY,largeL,'facecolor',[.5 .5 .5],'edgecolor','none'); light; colormap(gray); view([-20,30]); title('The set of tent poles') ```

Lower Bound Constraint Surface

The supporting poles determine a lower bound for the tent, L. We can visualize the constraint by plotting it as a magenta mesh.

```L = zeros(30); E = ones(2); L(15:16,15:16) = .5*E; L(5:6,5:6) = .3*E; L(25:26,5:6) = .3*E; L(5:6,25:26) = .3*E; L(25:26,25:26) = .3*E; % Add L to the plot. surface(L,'facecolor','none','edgecolor','m'); title('Lower Bound Constraint Surface') ```

Initial Guess for Optimization

To solve this problem, we will find the height of the optimized surface at a finite number of points. Our initial guess, `sstart`, is shown in blue.

```sstart = .5*ones(30,30); % Add it to the plot. surface(sstart,'FaceColor','none','LineStyle','none', ... 'Marker','.','MarkerEdgeColor','blue') title('Initial Value (blue) and Lower Bound (magenta)'); fig = gcf; fig.Renderer = 'zbuffer'; % Markers do not show up in OpenGL. ```

Resize Data for Optimization

In order to formulate the problem as a standard optimization problem, we resize both the matrices into vectors. `L` represents the initial values and `sstart` represents the lower boundary constraint.

```low = reshape(L,900,1); xstart = reshape(sstart,900,1); % Illustrate the reordering. % Draw grid points. xx = 0:4; [X, Y] = meshgrid(xx,xx); gpts = plot(X(:),Y(:),'b.'); gpts.MarkerSize = 10; axis off axis([-2 12 -1.5 5.5]) hold on % Draw arrow. l(1) = line([7.5 6.5],[2 2.5]); l(2) = line([7.5 6.5],[2 1.5]); l(3) = line([7.5 5.5],[2 2]); set(l,'color','b'); % Draw vector. yy = 0.2*xx; zz = [-1.5+yy,yy,1.5+yy,3+yy,4.5+yy]; vect = plot(9*ones(25,1),zz,'b.'); vect.MarkerSize = 9; axis off; hold off; ```

Formulate Optimization Problem

The surface formed by the elastic membrane is determined by the bound constrained optimization problem

` min{ c'*x+0.5*x'*H*x : low <= x }`

where `c'*x + 0.5*x'*H*x` is the discrete approximation of the energy function. `H` and `c` are as follows:

```H = delsq(numgrid('S',30+2)); h = 1/(30-1); c = -h^2*ones(30^2,1); ```

Hessian Sparsity Structure

Each point of the energy function is only affected by its immediate neighbors. Consequently the Hessian matrix `H` is sparse and has a special structure.

Because `H` is sparse we can use a large-scale algorithm to solve the optimization problem.

```spy(H); title('Sparsity Structure of Hessian Matrix'); ```

Run Optimization Solver

The trust region reflective algorithm in `quadprog` solves bound constrained quadratic problems. We select this algorithm by setting an option with `optimoptions`. Then solve the problem by calling `quadprog`.

```options = optimoptions('quadprog','Algorithm','trust-region-reflective',... 'Display','off'); x = quadprog(H, c, [], [], [], [], low, [], xstart, options); ```

Reshape Solution Back to Original Size

Now obtain the surface solution by going back to the original ordering using `reshape`. Then plot this solution in blue mesh.

```S = reshape(x,30,30); % Plot the starting surface. subplot(1,2,1); surf(L,'facecolor',[.5 .5 .5]); surface(sstart,'edgecolor','b','facecolor','none'); title('Starting Surface') axis off axis tight; view([-20,30]); % Plot the solution surface. subplot(1,2,2); surf(L,'facecolor',[.5 .5 .5]); surface(S,'edgecolor','b','facecolor','none'); title('Solution Surface') axis off axis tight; view([-20,30]); ```

Plot Circus Tent

Let's now plot the solution that the optimization solver found.

```subplot(1,1,1) surf(L,'facecolor',[0 0 0]); hold on; surfl(S); hold off; axis tight; axis off; view([-20,30]); ```