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

Рассмотрите создание цирка палатки, чтобы покрыть квадратный лот. Палатка имеет пять полюсов, покрытых тяжелым, упругим материалом. Задача состоит в том, чтобы найти естественную форму палатки. Моделируйте форму как высоту 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;

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

Решите проблему. Проигнорируйте полученный оператора «Ваш Гессиан не симметричен». 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]);