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

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')

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')

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.

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;

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

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');
```

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);

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]);

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]);

Was this topic helpful?