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

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

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

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;

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

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 функция. 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]);