Документация

Проблемное определение

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