Описание задачи

Полагайте, что создание цирка шапито покрывает квадратную партию. Палатка имеет пять полюсов, покрытых тяжелым, эластичным материалом. Проблема состоит в том, чтобы найти естественную форму палатки. Смоделируйте форму как высоту x (p) палатки в положении p.

Потенциальная энергия тяжелого материала, снятого к высоте x, является cx для постоянного c, который пропорционален весу материала. Для этой проблемы выберите c = 1/3000.

Эластичная потенциальная энергия части материала $$E_{stretch}$$ приблизительно пропорционально второй производной высоты материала, времена высота. Можно аппроксимировать вторую производную приближением конечной разности с пятью точками (примите, что шаги конечной разности имеют размер 1). Пусть $$\Delta x$$ представляйте сдвиг 1 в первом координатном направлении, и $$\Delta y$$ представляйте сдвиг 1 во втором координатном направлении.

Естественная форма палатки минимизирует общую потенциальную энергию. Путем дискретизации проблемы вы находите, что общая потенциальная энергия, чтобы минимизировать является суммой по всем положениям p $$E_{stretch}(p)$$ + cx (p).

Эта потенциальная энергия является квадратичным выражением в переменной x.

Задайте граничное условие, что высота палатки в ребрах является нулем. Крепления для палатки имеют сечение модуля 1 на 1, и палатка имеет общий размер 33 33 модулей. Задайте высоту и местоположение каждого полюса. Постройте квадратную область партии и крепления для палатки.

height = zeros(33);
height(6:7,6:7) = 0.3;
height(26:27,26:27) = 0.3;
height(6:7,26:27) = 0.3;
height(26:27,6:7) = 0.3;
height(16:17,16:17) = 0.5;
colormap(gray);
surfl(height)
axis tight
view([-20,30]);
title('Tent Poles and Region to Cover')

Сформулируйте задачу оптимизации

Создайте переменную x оптимизации представление высоты материала.

x = optimvar('x',size(height));

Установите x обнулять на контурах квадратной области.

boundary = false(size(height));
boundary([1,33],:) = true;
boundary(:,[1,33]) = true;
x.LowerBound(boundary) = 0;
x.UpperBound(boundary) = 0;

Вычислите эластичную потенциальную энергию в каждой точке. Во-первых, вычислите потенциальную энергию во внутренней части области, где конечные разности не переступают через область, содержащую решение.

L = size(height,1);
peStretch = optimexpr(L,L); % This initializes peStretch to zeros(L,L)
interior = 2:(L-1);
peStretch(interior,interior) = (-1*(x(interior - 1,interior) + x(interior + 1,interior) ...
    + x(interior,interior - 1) + x(interior,interior + 1)) + 4*x(interior,interior))...
    .*x(interior, interior);

Поскольку решение ограничивается быть 0 в ребрах области, вы не должны включать остаток от терминов. Все термины имеют кратное x, и x в ребре нуль. Для ссылки в случае, если вы хотите использовать различное граничное условие, следующее является прокомментированной версией потенциальной энергии.

% peStretch(1,interior) = (-1*(x(1,interior - 1) + x(1,interior + 1) + x(2,interior))...
%     + 4*x(1,interior)).*x(1,interior);
% peStretch(L,interior) = (-1*(x(L,interior - 1) + x(L,interior + 1) + x(L-1,interior))...
%     + 4*x(L,interior)).*x(L,interior);
% peStretch(interior,1) = (-1*(x(interior - 1,1) + x(interior + 1,1) + x(interior,2))...
%     + 4*x(interior,1)).*x(interior,1);
% peStretch(interior,L) = (-1*(x(interior - 1,L) + x(interior + 1,L) + x(interior,L-1))...
%     + 4*x(interior,L)).*x(interior,L);
% peStretch(1,1) = (-1*(x(2,1) + x(1,2)) + 4*x(1,1)).*x(1,1);
% peStretch(1,L) = (-1*(x(2,L) + x(1,L-1)) + 4*x(1,L)).*x(1,L);
% peStretch(L,1) = (-1*(x(L,2) + x(L-1,1)) + 4*x(L,1)).*x(L,1);
% peStretch(L,L) = (-1*(x(L-1,L) + x(L,L-1)) + 4*x(L,L)).*x(L,L);

Задайте потенциальную энергию из-за высоты материала, которая является x/3000.

peHeight = x/3000;

Создайте задачу оптимизации под названием tentproblem. Включайте выражение для целевой функции, которая является суммой двух потенциальных энергий по всем местоположениям.

tentproblem = optimproblem('Objective',sum(sum(peStretch + peHeight)));

Установите ограничение

Установите ограничение, что решение должно лечь выше значений height матрица. Эта матрица является нулем в большинстве местоположений, представляя землю, и включает высоту каждого крепления для палатки в его местоположении.

htcons = x >= height;
tentproblem.Constraints.htcons = htcons;

Запустите решатель оптимизации

Решите задачу. Проигнорируйте получившийся оператор "Your Hessian is not symmetric". solve вводит это выражение, потому что внутреннее преобразование от проблемной формы до квадратичной матрицы не гарантирует, что матрица симметрична.

sol = solve(tentproblem);

Solving problem using quadprog.
Your Hessian is not symmetric. Resetting H=(H+H')/2.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

Постройте решение

Постройте решение, найденное решателем оптимизации.

surfl(sol.x);
axis tight;
view([-20,30]);