В этом примере показано, как определить форму цирка шапито путем решения квадратичной задачи оптимизации. Палатка формируется из тяжелого, эластичного материала и приспосабливается к форме, которая имеет минимальную потенциальную энергию, удовлетворяющую ограничениям. Дискретизация проблемы приводит к связано ограниченной проблеме квадратичного программирования.
Для основанной на проблеме версии этого примера смотрите Связано ограниченное Квадратичное программирование, Основанное на проблеме.
Полагайте, что создание цирка шапито покрывает квадратную партию. Палатка имеет пять полюсов, покрытых тяжелым, эластичным материалом. Проблема состоит в том, чтобы найти естественную форму палатки. Смоделируйте форму как высоту x (p) палатки в положении p.
Потенциальная энергия тяжелого материала, снятого к высоте x, является cx для постоянного c, который пропорционален весу материала. Для этой проблемы выберите c = 1/3000.
Эластичная потенциальная энергия части материала приблизительно пропорционально второй производной высоты материала, времена высота. Можно аппроксимировать вторую производную приближением конечной разности с 5 точками (примите, что шаги конечной разности имеют размер 1). Пусть представляйте сдвиг 1 в первом координатном направлении, и представляйте сдвиг 1 во втором координатном направлении.
Естественная форма палатки минимизирует общую потенциальную энергию. Путем дискретизации проблемы вы находите, что общая потенциальная энергия, чтобы минимизировать является суммой по всем положениям 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
формулировка задачи должна минимизировать
.
В этом случае, линейный член соответствует потенциальной энергии высоты материала. Поэтому задайте f = 1/3000 для каждого компонента x.
f = ones(size(height))/3000;
Создайте матричное представление конечной разности при помощи delsq
функция. delsq
функция возвращает разреженную матрицу с записями 4 и-1 соответствия записям 4 и-1 в формуле для . Умножьте возвращенную матрицу на 2, чтобы иметь quadprog
решите квадратичную программу с энергетической функцией, как дано .
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]);