Связано ограниченное квадратичное программирование, основанное на решателе

В этом примере показано, как определить форму цирка шапито путем решения квадратичной задачи оптимизации. Палатка формируется из тяжелого, эластичного материала и приспосабливается к форме, которая имеет минимальную потенциальную энергию, удовлетворяющую ограничениям. Дискретизация проблемы приводит к связано ограниченной проблеме квадратичного программирования.

Для основанной на проблеме версии этого примера смотрите Связано ограниченное Квадратичное программирование, Основанное на проблеме.

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

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

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

Эластичная потенциальная энергия части материала Estretch приблизительно пропорционально второй производной высоты материала, времена высота. Можно аппроксимировать вторую производную приближением конечной разности с 5 точками (примите, что шаги конечной разности имеют размер 1). Пусть Δx представляйте сдвиг 1 в первом координатном направлении, и Δy представляйте сдвиг 1 во втором координатном направлении.

Estretch(p)=(-1(x(p+Δx)+x(p-Δx)+x(p+Δy)+x(p-Δy))+4x(p))x(p).

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

Figure contains an axes. The axes with title Tent Poles and Region to Cover contains an object of type surface.

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

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 формулировка задачи должна минимизировать

12xTHx+fTx.

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

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

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

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

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

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

Figure contains an axes. The axes with title Sparsity Structure of Hessian Matrix contains an object of type line.

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

Решите задачу путем вызова 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]);

Figure contains an axes. The axes contains an object of type surface.

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте