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

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

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

Упругая потенциальная энергия части материала $$E_{stretch}$$ приблизительно пропорционально второй производной от высоты материала, умноженной на высоту. Можно аппроксимировать вторую производную 5-точечного конечного различия приближения (предположим, что конечные шаги различия имеют размер 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')

Создание граничных условий

The height матрица определяет нижние границы решения x. Чтобы ограничить решение нулем на контуре, установите верхнюю границу ub чтобы быть нулем на контуре.

boundary = false(size(height));
boundary([1,33],:) = true;
boundary(:,[1,33]) = true;
ub = inf(size(boundary)); % No upper bound on most of the region
ub(boundary) = 0;

Создайте матрицы целевых функций

The quadprog формулировка задачи состоит в том, чтобы минимизировать

$$\frac{1}{2}{x}^{T}Hx+f^{T}x$$.

В этом случае линейный термин $$f^{T}x$$ соответствует потенциальной энергии высоты материала. Поэтому задайте f = 1/3000 для каждого компонента x.

f = ones(size(height))/3000;

Создайте конечную матрицу различий, представляющую $$E_{stretch}$$ при помощи delsq функция. The delsq функция возвращает разреженную матрицу с записями 4 и -1, соответствующими записям 4 и -1 в формуле для $$E_{stretch}(p)$$. Умножьте возвращенную матрицу на 2, чтобы иметь quadprog решить квадратичную программу с энергетической функцией, заданной $$E_{stretch}$$.

H = delsq(numgrid('S',33+2))*2;

Просмотрите структуру матрицы H. Матрица работает с x(:), что означает матрицу x преобразуется в вектор линейным индексированием.

spy(H);
title('Sparsity Structure of Hessian Matrix');

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

Решите проблему, позвонив quadprog.

x = quadprog(H,f,[],[],[],[],height,ub);

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.

Решение для построения графика

Измените форму решения x в матрицу S. Затем постройте график решения.

S = reshape(x,size(height));
surfl(S);
axis tight;
view([-20,30]);