Определение проблемы

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

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

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

Естественная форма палатки минимизирует общую потенциальную энергию. Дискретизируя задачу, вы обнаруживаете, что общая потенциальная энергия для минимизации является суммой по всем позициям p $${}_{\mathrm{Estretch}}(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]);