Определение проблемы

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

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

Энергия упругого потенциала куска материала $${}_{\mathrm{Estretch}}$$ приблизительно пропорциональна второй производной высоты материала, умноженной на высоту. Вторую производную можно аппроксимировать аппроксимацией 5-точечной конечной разности (предположим, что шаги конечной разности имеют размер 1). Пусть $$\mathrm{\Delta x}$$ представляет сдвиг 1 в первом направлении координат, а $$\mathrm{\Delta y}$$ представляет сдвиг на 1 во втором направлении координат.

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

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

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{\mathrm{}}{\mathrm{}}{}^{}\mathrm{12xTHx}{+}^{}$$fTx.

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

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

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

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]);