www.gusucode.com > optim 案例源码 matlab代码程序 > optim/circustent.m
%% 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(TM). % Copyright 1990-2014 The MathWorks, Inc. %% 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]);