Метод градиентного спуска

Производная функционала

\(Y\) — пространство состояний, \(U\) — пространство управлений, \(V\) — пространство правых частей, \(U_{\mathrm{ad}} \subset U\) — множество допустимых управлений.

\(F \colon Y \times U \to V\) — оператор ограничений, \(J \colon Y \times U \to \mathbb R\) — функционал качества.

Постановка задачи: \(J(y,u) \to \inf,\quad F(y,u) = 0, \;\; u \in U_{\mathrm{ad}}.\)

Также рассмотрим приведенную задачу. Составим приведенный функционал качества: \(\widehat J(u) = J(y(u), u),\) где \(y(u)\) — решение уравнения \(F(y,u) = 0\). Приведенная задача: \(\widehat J(u) \to \inf,\quad u \in U_{\mathrm{ad}}.\)

Для применения градиентного метода необходимо вычислить производную Фреше приведенного функционала \(\widehat J'(u) \in U'\). Если \(U\) — гильбертово пространство, то \(U'\) можно отождествить с \(U\) и считать, что \(\widehat J'(u) \in U\).

Составим функцию Лагранжа \(L \colon Y \times U \times V' \to \mathbb R\): \(L(y,u,p) = J(y,u) + \langle p, F(y,u) \rangle_{V' \times V}.\)

Подставим в функцию Лагранжа \(y(u)\) вместо \(y\): \(L(y(u), u, p) = J(y(u), u) + \langle p, F(y(u), u) \rangle_{V'\times V} = J(y(u), u) = \widehat J(u),\) поэтому \(\widehat J(u) = L(y(u), u, p).\) В этом равенстве \(p \in V'\) — произвольный элемент.

Теперь продифференцируем это равенство по \(u\): \(\widehat J'(u) = L_y(y(u), u, p)y'(u) + L_u(y(u), u, p).\)

Выберем специальным образом \(p = p(u) \in V'\) такое, что \(L_y(y(u), u, p) = 0.\) (Это сопряженное уравнение.) Тогда \(\widehat J'(u) = L_u(y(u), u, p(u)).\)

Алгоритм

Если \(U_{\mathrm{ad}} = U\), то метод градиентного спуска принимает вид: \(u^{k+1} = u^k - \lambda \widehat J'(u^k),\) где \(\lambda > 0\) — итерационный параметр. \(u^0\) — некоторое начальное приближение.

Если же \(U_{\mathrm{ad}} \neq U\), то применяем метод проекции градиента: \(u^{k+1} = P_{U_{\mathrm{ad}}}(u^k - \lambda \widehat J'(u^k)),\) где \(P_{U_{\mathrm{ad}}}(u)\) — проекция \(u \in U\) на множество \(U_{\mathrm{ad}}\).

Пример

\(\tag{pr1} y_t = y_{xx} + u,\quad x \in (0,1),\;\; t \in (0,T),\) \(\tag{pr2} y|_{x=0;1} = 0,\) \(\tag{pr3} y|_{t=0} = y_0(x),\quad x \in \Omega.\) \(J(y,u) = \frac{1}{2}\int_\Omega(y|_{t=T} - y_d)^2dx + \frac{N}{2}\int_0^T \int_\Omega u^2 dxdt \to \inf.\)

Вычислим функцию Лагранжа (подробности см. в лекциях): \(L(y,u,p,q) = J(y,u) + \int_0^T(y_t + Ay - u, p)dt + (q, y|_{t=0} - y_0).\)

Возьмем производную по \(u\): \(\langle L_u(y,u,p,q), v \rangle_{U' \times U} = N\int_0^T(u, v)dt - \int_0^T (p, v)dt = \int_0^T(Nu - p, v)dt.\) Поэтому \(\widehat J'(u) = L_u(y,u,p(u),q(u)) = Nu - p(u).\) Здесь \(p(u)\) — решение сопряженного уравнения при заданном \(u\). Таким образом, для вычисления градиента приведеного функционала \(\widehat J'(u)\) необходимо решить сопряженную задачу при заданном \(u\).

Сопряженная система: \(\tag{sopr1} -p_t - p_{xx} = 0,\) \(\tag{sopr2} p|_{x=0;1} = 0,\) \(\tag{sopr3} p|_{t=T} = y_d - y|_{t=T}.\)

Таким образом, для вычисления следующего приближения \(u^{k+1}\) нужно сначала решить прямую задачу (pr1)–(pr3), затем подставить найденное решение \(y\) при \(t=T\) в сопряженную задачу (sopr1)–(sopr3) и решить ее. После этого следующее приближение находится по формуле: \(u^{k+1} = u^k - \lambda(Nu^k - p^k) = (1 - \lambda N)u^k + \lambda p^k.\) Здесь \(p^k = p(u^k)\). Заметим, что при \(\lambda = 1/N\) градиентный метод переходит в метод простой итерации. Чтобы метод градиентного спуска сходился, нужно выбрать малое значение \(\lambda\), но не слишком малое, чтобы скорость сходимости была не слишком медленной.

Программная реализация

Вводим сетку: \(x_i = ih,\;\; i=\overline{0,M},\;\; h = 1/M,\) \(t_k = k\tau,\;\;k=\overline{0,K},\;\; \tau = T/K.\)

Чисто неявная схема: \(\tag{rs1} \frac{y_i^{k+1} - y_i^k}{\tau} = \frac{y_{i+1}^{k+1} - 2y_i^{k+1} + y_{i-1}^{k+1}}{h^2} + u_i^{k+1}, \quad i = \overline{1,M-1},\;\; k = \overline{0,K-1},\) \(y_0^{k+1} = y_M^{k+1} = 0,\) \(y_i^0 = y_0(x_i).\)

Сопряженная задача сводится к обычному уравнению теплопроводности заменой \(\widetilde p(t) = p(T-t)\): \(\widetilde p_t - \widetilde p_{xx} = 0,\) \(\widetilde p|_{x=0;1} = 0,\) \(\widetilde p|_{t=0} = y_d - y|_{t=T}.\)

Разностная схема: \(\tag{rs2} \frac{p_i^{k} - p_i^{k+1}}{\tau} = \frac{p_{i+1}^k - 2p_i^k + p_{i-1}^k}{h^2}, \quad i = \overline{1,M-1},\;\; k = \overline{0, K-1},\) \(p_0^k = p_M^k = 0,\) \(p_i^K = y_d(x_i) - y_i^K.\)

Необходимо сделать несколько итераций градиентного спуска. Условием окончания итераций может служить \(\|u^{k+1} - u^k\| < \varepsilon.\)

Запишем (rs1) в виде: \(-\frac{1}{h^2}y_{i-1}^{k+1} + \left(\frac{1}{\tau} + \frac{2}{h^2} \right)y_i^{k+1} - \frac{1}{h^2}y_{i+1}^{k+1} = u_i^{k+1} + \frac{y_i^k}{\tau}.\)

Аналогично запишем (rs2) в виде: \(-\frac{1}{h^2}p_{i-1}^{k} + \left(\frac{1}{\tau} + \frac{2}{h^2} \right)p_i^{k} - \frac{1}{h^2}p_{i+1}^{k} = \frac{p_i^{k+1}}{\tau}.\)

Исходя из этих формул, можно составить матрицу и вектор правой части для СЛАУ \(Ay = f\). В обоих случаях матрица будет одинаковой: \(a_{ii} = \frac{1}{\tau} + \frac{2}{h^2},\quad a_{i,i-1} = a_{i,i+1} = -\frac{1}{h^2},\quad i = \overline{1,M-1}.\) (Остальные элементы равны 0.)

Правые части: \(b_i = u_i^{k+1} + \frac{y_i^k}{\tau}, \quad i = \overline{1,M-1},\) \(c_i = \frac{p_i^{k+1}}{\tau}, \quad i = \overline{1,M-1}.\)

Но также нужно заполнить 0-ю и \(M\)-ю строки. \(a_{ii} = 1, \quad b_i = c_i = 0, \quad i = 0,M.\)

Для реализации можно использовать свободно распространяемую среду Octave или облако программирования Wolfram: http://www.wolfram.com/programming-cloud.

В Octave индексы начинаются с 1, а не с 0, поэтому нужно будет прибавить единицу к индексам в этих формулах.

Информацию о среде Octave можно получить на сайтах: http://gnu-octave.narod.ru, http://elis.dvo.ru/~shamray/manual/gnu.pdf, http://www.altlinux.org/images/0/07/OctaveBook.pdf, а также в книге Hansen J.S. GNU Octave. Beginner’s Guide. 2011.

Скачать Octave можно по адресу: http://goo.gl/2KOI7c.

Если вы планируете использовать дополнительные пакеты, необходимо поставить в установщике галочку “Octave Forge”.

Чтобы запустить программу, создайте файл с расширением .m (например, opt.m), наберите команду cd путь_к_папке и вызовите команду opt (или ваше имя файла).

Матрицу \(A\) объявляем как разреженную (ключевое слово sparse), чтобы Octave использовала эффективные алгоритмы решения СЛАУ с разреженными матрицами (например, метод прогонки). Пример см. в книге Хансена, с. 180.

T = 10; % временной промежуток
N = 0.01; % параметр регуляризации

lambda = 20; % итерационный параметр
numsteps = 50; % количество итераций

function yd = yd(x) % желаемая температура при t = T
  yd = 0.5;
endfunction

function y0 = y0(x) % начальная температура
  y0 = sin(pi*x);
endfunction

% параметры сетки
M = 100;
K = 100;
global h = 1/M;
tau = T/K;

function xi = xi(i)
  global h;
  xi = (i-1)*h;
endfunction

u = zeros(M+1, K+1); % инициализация управления

% заполняем матрицу СЛАУ
A = sparse(M+1, M+1);
A(1, 1) = A(M+1, M+1) = 1;
for i = 2:M
  A(i, i-1) = A(i, i+1) = -1/h^2;
  A(i, i) = 1/tau + 2/h^2;
endfor

for step = 1:numsteps

  % начальные условия
  for i = 1:M+1
    y(i, 1) = y0(xi(i));
  endfor

  %решаем прямую задачу
  b = zeros(M+1, 1);
  for k = 1:K
    % заполняем вектор правой части
    b(1) = b(M+1) = 0;
    for i = 2:M
      b(i) = u(i, k+1) + y(i, k)/tau;
    endfor
    y(:, k+1) = A\b;
  endfor

  % температура при t = T
  %figure
  %x = linspace(0, 1, M+1);
  %plot(x, y(:, K+1));

  % температура при всех t
  %figure
  %x = linspace(0, 1, M+1);
  %t = linspace(0, T, K+1);
  %[TT XX] = meshgrid(t, x);
  %surf(TT, XX, y);
  %view(50, 30); % угол обзора

  %терминальные условия
  for i = 1:M+1
    p(i, K+1) = yd(xi(i)) - y(i, K+1);
  endfor

  %решаем сопряженную задачу
  c = zeros(M+1, 1);
  for k = K:-1:1
    % заполняем вектор правой части
    c(1) = c(M+1) = 0;
    for i = 2:M
      c(i) = p(i, k+1)/tau;
    endfor
    p(:, k) = A\c;
  endfor

  % сопряженная температура при всех t
  %figure
  %x = linspace(0, 1, M+1);
  %t = linspace(0, T, K+1);
  %[TT XX] = meshgrid(t, x);
  %surf(TT, XX, p);
  %view(50, 30); % угол обзора

  % пересчитываем управление
  %for k = 1:K+1
  % for i = 1:M+1
  % u1(i, k) = (1 - lambda*N)*u(i, k) + lambda*p(i, k);
  % endfor
  %endfor
  u1 = (1 - lambda*N)*u + lambda*p;

  % вычисляем норму разности приближений
  d = 0;
  for k = 1:K+1
    for i = 1:M+1
      d += (u(i, k) - u1(i, k))^2;
    endfor
  endfor
  d /= (K+1)*(M+1);
  d = sqrt(d);

  u = u1;

endfor % step

d % выводим разность норм
% управление
figure
x = linspace(0, 1, M+1);
t = linspace(0, T, K+1);
[TT XX] = meshgrid(t, x);
surf(TT, XX, u);
view(50, 30); % угол обзора

% температура при t = T
figure
x = linspace(0, 1, M+1);
y1 = arrayfun(@yd, x);
plot(x, y(:, K+1), x, y1);

Метод Фурье

Метод Фурье применим, вероятно, только для вариантов, в которых нет дополнительного ограничения на управление (когда, например, управление равно 0 на части интервала).

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