Оптимальное управление — 07 — Конспект

$\global\def\at#1{\left. #1 \right\rvert}$ $\global\def\abs#1{\left\lvert #1 \right\rvert}$ $\global\def\norm#1{\left\lVert #1 \right\rVert}$ $\global\def\bvec#1{\mathbf{#1}}$ $\global\def\floor#1{\left\lfloor #1 \right\rfloor}$ $\global\def\limto#1{\underset{#1}{\longrightarrow}}$ $\global\def\prob#1{\mathbb{P} \left\{ #1 \right\}}$ $\global\def\mean#1{\mathbb{E} \left[ #1 \right]}$ $\global\def\disp#1{D \left[ #1 \right]}$ $\global\def\dp#1#2{#1 \cdot #2\,}$ $\global\def\vp#1#2{#1 \times #2\,}$ $\global\def\dv#1#2{\frac{d #1}{d #2}}$ $\global\def\rdv#1#2{\frac{d' #1}{d #2}}$ $\global\def\pd#1#2{\frac{\partial #1}{\partial #2}}$ $\global\def\pdv2#1#2{\frac{\partial^2 #1}{\partial #2^2}}$ $\global\def\pdvk#1#2#3{\frac{\partial^#1 #2}{\partial #3^#1}}$ $\global\def\ppdv#1#2#3{\frac{\partial^2 #1}{\partial #2 \partial #3}}$ $\global\def\pois#1{\left\{ #1 \right\}}$ $\global\def\paren#1{\left( #1 \right)}$ $\global\def\bydef#1{\overset{\mathrm{def}}{#1}}$ $\global\def\mbox#1{\text{#1}}$ $\global\def\div{\text{div}\,}$ $\global\def\dsum{\displaystyle\sum\,}$ $\global\def\grad{\text{grad}\,}$ $\global\def\rot{\text{rot}\,}$ $\global\def\vb#1{\textbf{#1}}$ $\global\def\op#1{\mathrm{#1}\,}$ $\global\def\Im{\text{Im}\,}$ $\global\def\Res{\text{Res}\,}$ $\global\def\Re{\text{Re}\,}$ $\global\def\argtg{\text{argtg}\,}$ $\global\def\ch{\text{ch}\,}$ $\global\def\const{\text{const}\,}$ $\global\def\degree{\text{degree}\,}$ $\global\def\proj{\mathrm{proj}}$ $\global\def\rank{\mathrm{rank}}$ $\global\def\res{\text{res}\,}$ $\global\def\sh{\text{sh}\,}$ $\global\def\sign{\text{sign}\,}$ $\global\def\tg{\mathrm{tg}\,}$

2023-09-14

Глава 1. Введение в ПМ

Постановка задачи

  1. уравнение движения: \[ \dot x(t) = f(x(t), u(t), t), \] где $x \in \mathbb{R}^n, u \in \mathbb{R}^m, t \in \mathbb{R}$;
  2. критерий качества:
    1. Лагранжа: \[ J_0 = \int\limits_{t_\text{н}}^{t_\text{к}} F_0(x(t), u(t), t) dt, \qquad F_0 \in \mathbb{R} \]
    2. Майера: \[ J_0 = \varphi_0(x(t_\text{к})), \qquad \varphi_0 \in \mathbb{R} \]
    3. Больца \[ J_0 = \int\limits_{t_\text{н}}^{t_\text{к}} F_0(x(t), u(t), t) dt + \varphi_0(x(t_\text{к})) \]
  3. (опц.) ограничения
    1. траектория
      • изопериметрическая \[ J_i = \int\limits_{t_\text{н}}^{t_\text{к}} F_i(x(t)) dt \leqslant 0 \] или $\geqslant, =$.
      • с фиксированным концом \[ x(t_\text{н}) = x_\text{н}, \qquad x(t_\text{к}) = x_\text{к}. \]
      • со свободным концом \[ x(t_\text{н}) = x_\text{н} \quad \text{или} \quad x(t_\text{к}) = x_\text{к}. \]
      • с плавающими концами \[ x(t_\text{н}) \in G(t_\text{н}), \qquad x(t_\text{к}) \in G(t_\text{к}). \]
    2. управление
      • программное управление \[ u = u(t) \]
        ПМП можно применять только для этого типа управления.
      • $u \in U(t)$ — область может меняться со временем
      • изопериметрическое ограничение: \[ J_j = \int\limits_{t_\text{н}}^{t_\text{к}} F_j(u(t), t) dt = 0, \] либо $\geqslant, \leqslant$.
(оптимального управления).
Формально можно все функционалы представить в виде Больца: \[ J_k = \int\limits_{t_\text{н}}^{t_\text{к}} F_k(x(t), u, t) dt + \varphi_k(t_\text{н}, t_\text{к}, x(t_\text{н}), x(t_\text{к})), \quad k = \overline{0, m_2}, \] где \[ \begin{gathered} F_k \in \mathbb{R}, \quad \varphi_k \in \mathbb{R}, \\ u(t) \in PC, \quad x(t) \in PCC, \end{gathered} \] причём $PC$ — кусочно-непрерывные функции, $PCC$ — кусочно-непрерывно дифференцируемые функции.
Рассмотрим совокупность варьируемых параметров: \[ q = q(t_\text{н}, t_\text{к}, x(t_\text{н}), x(t_\text{к}), u(t)). \]
Совокупность $q$ варьируемых параметров называется допустимой, если она удовлетворяет условиям \[ \begin{aligned} \dot x(t) &= f(x(t), u(t), t), \\ J_k &= \int\limits_{t_\text{н}}^{t_\text{к}} F_k(x(t), u, t) dt + \varphi_k(t_\text{н}, t_\text{к}, x(t_\text{н}), x(t_\text{к})). \end{aligned} \]
Совокупность $\tilde q$ варьируемых параметров называется оптимальной, если для любой допустимой совокупности $q$, близкой к $\tilde q$, выполняется соотношение \[ J_0(\tilde q) \leqslant J_0(q). \]

2023-09-21

Проварьируем функцию $x(t)$: \[ \tilde x(t) = x(t) + \varepsilon \xi(t), \qquad \xi(t) \in C^1. \] Тогда изохронной вариацией будем называть \[ \delta x \bydef= \tilde x(t) - x(t) = \varepsilon \xi(t). \] При $\varepsilon \to 0$ вариация также будет стремиться к нулю: $\delta x \to 0$, то есть будет бесконечно малой.

Вариацию \[ \Delta x \bydef= \tilde x(t + \Delta t) - x(t) \] будем называть бесконечно малой полной вариацией. Связь полной и изохронной вариациями выражается так: \[ \Delta x = \delta x + \dot x \Delta t. \]

Пусть $\Phi(x(t), t)$ — функционал, тогда его полная вариация представляется в виде \[ \Delta \Phi(x(t), t) = \Phi(\tilde x(t + \Delta t), t + \Delta t) - \Phi(x(t), t). \] Как и для функции $x(t)$, для функционала $\Phi(x(t), t)$ справедлива следующая связь: \[ \Delta \Phi = \delta \Phi + \dot \Phi \Delta t. \]

Рассмотрим задачу оптимального управления:

Будем считать, что целевой функционал представлен в виде Больца \[ J_0 = \int\limits_{t_\text{н}}^{t_\text{к}} F_0 dt + \varphi_0. \]

Предполагаем, что

  1. $x(t)$ — (кусочно-) непрерывно-дифференцируемая функция;
  2. $u(t)$ — (кусочно-) непрерывная функция;
  3. $\displaystyle \pd{F_0}{x} \in C^1$.

Тогда, введя множители Лагранжа:

Поставленную задачу можно свести к задаче безусловной минимизации, составив условный функционал: \[ \Psi = \int\limits_{t_\text{н}}^{t_\text{к}} \left( F_0 + \lambda(\dot x - f) \right) dt + \varphi_0 + \sum\limits_{j=1}^{l} \nu_j J_j. \] Введём обозначения: \[ L = F_0 + \lambda (\dot x - f), \qquad R = \varphi_0 + \sum\limits_{j=1}^{l} \nu_j J_j. \] Тогда \[ \Psi = \int\limits_{t_\text{н}}^{t_\text{к}} L dt + R \to \min. \]

Необходимое условие оптимальности: \[ \Delta \Psi = 0. \]

Введём в рассмотрение функцию Гамильтона: \[ H = H(t, x, \lambda, u) = -F_0 + \lambda f. \]

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

Если $u(t)$ — непрерывная функция без ограничений, тогда следующие условия являются необходимыми условиями оптимальности:

  1. уравнения движения: \[ \dot x_i = \pd{H}{\lambda_i} = f_i(x, u, t); \]
  2. уравнения Эйлера-Лагранжа: \[ \dot \lambda_i = - \pd{H}{x_i} = \pd{F_0}{x_i} - \sum\limits_{j=1}^{n} \lambda_j \pd{f_j}{x_i}; \]
  3. условия трансверсальности: \[ \Delta R(t_\text{н}, t_\text{к}, x(t_\text{н}), x(t_\text{к})) - \left. \paren{\lambda \Delta x - H \Delta t} \right|_{t_\text{н}}^{t_\text{к}} = 0; \]
  4. $\displaystyle \pd{H}{u} = 0$.
(необходимые условия оптимальности + ПМП).

Если $u$ — кусочно-непрерывная функция с заданными ограничениями ($u \in U$), то необходимыми условиями оптимальности являются следующие условия:

  1. уравнения движения: \[ \dot x_i = \pd{H}{\lambda_i} = f_i(x, u, t); \]
  2. уравнения Эйлера-Лагранжа: \[ \dot \lambda_i = - \pd{H}{x_i} = \pd{F_0}{x_i} - \sum\limits_{j=1}^{n} \lambda_j \pd{f_j}{x_i}; \]
  3. условия трансверсальности: \[ \Delta R(t_\text{н}, t_\text{к}, x(t_\text{н}), x(t_\text{к})) - \left. \lambda \Delta x - H \Delta t \right|_{t_\text{н}}^{t_\text{к}} = 0; \]
  4. $\displaystyle H(t, x, \lambda, u) = \max_{u' \in U} H(t, x, \lambda, u')$.

Сформулируем ПМП.

Рассмотрим уравнения движения \[ \dot x(t) = f(x(t), u(t), t), \qquad x \in \R^n, \quad u \in \R^m, \quad t \in [t_\text{н}, t_\text{к}] \] с начальными условиями $x(t_\text{н}) = x^\text{н}$.

В качестве критерия качества рассмотрим задачу Больца со свободным правым концом \[ J_0 = \int\limits_{t_\text{н}}^{t_\text{к}} F_0(x(t), u(t), t) dt + \varphi_0(x(t_\text{к})) \to \inf. \]

Управление $u(t) \in U$ будем искать в классе кусочно-непрерывных функций, то есть $u(t) \in PC$.

(принцип максимума Понтрягина).

Если $u_0(t)$ — оптимальное управление для поставленной задачи, а $x_0(t)$ — соответствующая ему оптимальная траектория, то $u_0$ удовлетворяет условиям максимальности: \[ \max_{u \in U} H(t, x_0(t), \lambda(t), u(t)) = H(t, x_0(t), \lambda(t), u_0(t)), \] где \[ H(t, x(t), \lambda(t), u(t)) = -F_0(x(t), u(t), t) + \lambda(t) f(x(t), u(t), t), \] а множители $\lambda(t) \in C, \lambda(t) \neq 0$ являются решением уравнений Эйлера-Лагранжа \[ \dot \lambda = - \pd{H}{x} \] и удовлетворяют условиям трансверсальности: \[ \lambda(t_\text{к}) = - \pd{\varphi_0(x_0(t_\text{к}))}{x_\text{к}}. \]

В общем случае ПМП является необходимым условием оптимальности, однако он является достаточным для некоторых классов задач (например, для линейно выпуклых — где $f$ линейна по $x$, а $J_0$ — выпуклый)

2023-09-28

Посмотреть задачи в группе в тимсе.

Глава 2. Оптимальное управление механическими системами

Оптимальное по расходу топлива гашение плоских колебаний спутника

  1. Постановка задачи.
    Введём системы координат: условно неподвижная $C xyz$ и связанная со спутником $C x_0 y_0 z_0$.
    Движение $C xyz$ относительно $C x_0 y_0 z_0$ определяется углами Эйлера. Колебания плоские, поэтому $Cy = Cy_0$. Значит, в плоскости орбиты есть угол отклонения $\theta$.

    Уравнение, описывающее вращение тела относительно неподвижной оси: \[ I_y \ddot \theta = \sum\limits_{i=1}^{N} M_y(F_i^{(e)}) = M_y. \] где $F_i^{(e)}$ — внешние силы (индекс сверху означает это), $I_y$ — момент инерции спутника относительно оси $y$, $M_y$ — момент внешней силы относительно оси $y$.

    Выразим проекцию $M_y$ момента внешних сил на $y$: \[ M_y = -3 \omega_0^2 (I_x - I_z) \sin\theta \cos\theta, \] где $\omega_0$ — скорость орбитального движения: \[ \omega_0 = \kappa r_c^{-3/2} = \const, \] где $\kappa = G m$, где $m$ — масса планеты, $G$ — гравитационная постоянная.

    Будем предполагать, что $I_x \gt I_z$ (то есть спутник вытянут относительно оси $z$), отсюда \[ I_y \ddot \theta + 3 \omega_0^2 (I_x - I_z) \sin\theta \cos\theta = 0, \] причём $(I_x - I_z) \gt 0$.

    Предполагаем также, что $\omega_0 \gt 0$. Заметим, момент инерции — величина положительная, поэтому введём обозначение \[ \omega^2 := 3 \omega_0^2 (I_x - I_z) I_y^{-1} \gt 0. \] Тогда \[ \ddot \theta + \omega^2 \sin\theta \cos\theta = 0. \] Пусть $\varphi = 2 \theta$, тогда \[ \ddot \varphi + \omega^2 \sin\varphi = 0. \] Отклонения малые, поэтому $\sin\varphi \sim \varphi$, и получаем уравнение математического маятника: \[ \ddot \varphi + \omega^2 \varphi = 0. \]

    Предполагаем, что топливо расходуется настолько мало, что моменты инерции можно считать постоянными, то есть $\omega^2 = \const$.

    Введём управление. Предполагаем, что на осях располагается пара двигателей, причём одна из пар создаёт положительный момент, другая — отрицательный момент. Обозначив положительный момент как $u_1 \gt 0$, отрицательный — как $u_2 \gt 0$, приходим к уравнению \[ \ddot \varphi + \omega^2 \varphi = u_1 - u_2. \] Моменты создаются микродвигателями (считаем, что они одинаковые), поэтому они ограничены: \[ 0 \leqslant u_i \leqslant \nu_i, \qquad i = 1,2, \quad \nu_1 = \nu_2 = \nu = \const \gt 0. \]

    Введём ограничения. Предполагаем, что в начальный момент $t_\text{н} = 0$ времени было отклонение, то есть \[ \varphi(0) = \psi_0, \qquad \dot \varphi(0) = \dot \psi_0. \] В конечный момент времени $\varphi(t_\text{к}) = 0, \quad \dot \varphi(t_\text{к}) = 0$.

    Введём критерий качества. Приведём задачу к каноническому виду (т.е. к системе ОДУ первого порядка): \[ x_1 = \varphi, \qquad x_2 = \dot \varphi. \] Тогда \[ \begin{equation} \tag{1} \begin{aligned} \dot x_1 &= x_2, \\ \dot x_2 &= - \omega^2 x_1 + u_1 - u_2, \end{aligned} \end{equation} \] причём \[ \begin{aligned} x_1(0) &= x_1^0, & x_2(0) &= x_2^0, \\ x_1(t_\text{к}) &= 0, & x_2(t_\text{к}) &= 0, \end{aligned} \tag{2} \] \[ 0 \leqslant u_{1,2} \leqslant \nu \tag{3} \]

    Условия трансверсальности работать не будут, тк ограничения на траекторию мы уже расписали.

    По теореме Калмана проверяется, что система управляема.

    Рассмотрим 2 критерия качества: расход топлива \[ \tag{4} J_1 = \int\limits_{t_\text{н}}^{t_\text{к}} (u_1 + u_2) dt \] и задача быстродействия \[ \tag{5} J_2 = \int\limits_{t_\text{н}}^{t_\text{к}} dt = t_\text{к} - t_\text{н}. \]

    Почему $J_1$ именно такого вида? Покажем, что именно такой вид пропорционален расходу топлива. Пусть реактивная тяга одного микродвигателя обозначается как $\mu u_r$, где $u_r$ — постоянная относительная (относительно спутника) скорость истечения топлива, а $\mu = - \dv{m}{t}$ — расход топлива.

    Пусть $\mu_1$ — расход топлива двух микродвигателей, создающих положительный момент. Понятно, что $\mu_1 = 2 \mu$. Тяга каждого микродвигателя тогда получается как $\dfrac{\mu_1}{2} u_r$. Пусть плечо $l$, тогда момент равен \[ m^+ = \dfrac{\mu_1}{2} u_r l. \] Считаем, что \[ u_1 = \frac{m^+}{I_y}, \qquad u_2 = \frac{m^-}{I_y}. \] Тогда \[ \tag{6} \begin{aligned} u_1 &= \frac{\mu_1}{2 I_y} u_r l, \\ u_2 &= \frac{\mu_2}{2 I_y} u_r l. \end{aligned} \] Подставив эти значения в критерий качества $J_1$, получим, что \[ \begin{aligned} J_1 &= \int\limits_{t_\text{н}}^{t_\text{к}} \paren{ \frac{\mu_1}{2 I_y} u_r l + \frac{\mu_2}{2 I_y} u_r l } dt = \\ &= \frac{u_r l}{2 I_y} \int\limits_{t_\text{н}}^{t_\text{к}} \paren{\mu_1 + \mu_2} dt. \end{aligned} \] Заметим, что $\mu = \mu_1 + \mu_2 = -\dv{m}{t}$ — полный расход топлива. Проинтегрируем, тогда \[ J_1 = \frac{u_r l}{2 I_y} (m_\text{н} - m_\text{к}). \]

  2. Исследование функции Гамильтона.
    Составим функцию Гамильтона: \[ \tag{7} \begin{aligned} H &= -F_0 + \lambda f = \\ &= - (u_1 + u_2) + \lambda_1 x_2 + \lambda_2 (-\omega^2 x_1 + u_1 - u_2). \end{aligned} \] Выпишем уравнения Эйлера-Лагранжа \[ \tag{8} \begin{aligned} \dot \lambda_1 = - \pd{H}{x_1} = \omega^2 \lambda_2 \\ \dot \lambda_2 = - \pd{H}{x_2} = - \lambda_1. \end{aligned} \] Тогда \[ \ddot \lambda_2 = - \dot \lambda_1 = - \omega^2 \lambda_2, \] то есть \[ \ddot \lambda_2 + \omega^2 \lambda_2 = 0. \] Выпишем решение: \[ \begin{aligned} \lambda_1 &= - a \omega \cos(\omega t + \alpha), \\ \lambda_2 &= \phantom{-} a \sin(\omega t + \alpha), \end{aligned} \] где $a = \const \gt 0$ — амплитуда, $\alpha = \const$ — начальная фаза колебания (находится из начальных условий).
    Так как $\lambda_1, \lambda_2$ имеют периодический характер, то и управление будет иметь такой же характер.

    Понятно, что амплитуда $a \neq 0$ (иначе $\lambda = 0$, что невозможно).

    Вырожденный режим \[ \pd{H}{u_1} = \pd{H}{u_2} = 0 \] в нашем случае невозможен. Следовательно, нет ни одного временного промежутка, где оба управления являются оптимальными, поэтому эти равенства определяют точки переключения соответствующего управления.

    \[ \begin{gathered} \pd{H}{u_1} = 0, \implies \lambda_2 = 1 \implies \\ \implies a \sin(\omega t + \alpha) = 1 \end{gathered} \] тогда \[ \omega t + \alpha = \frac{\pi}{2} \mp \arccos (1/a) \gt 0, \implies a \gt 1. \] Обозначим $\sigma = \arccos(1/a)$, тогда \[ \omega t + \alpha = \frac{\pi}{2} \mp \sigma + 2 \pi n. \] Это — точки переключения первого управления. Если $-\sigma$, то точка включения, если $\sigma$ — выключения: \[ \begin{aligned} \omega \tau_1 + \alpha &= \frac{\pi}{2} - \sigma + 2 \pi n, \\ \omega \tau_2 + \alpha &= \frac{\pi}{2} + \sigma + 2 \pi n, & n \in \mathbb{N}. \end{aligned} \] Можно записать как \[ \tag{9} \begin{aligned} \omega \tau_1 + \alpha &= \frac{\pi}{2} (4 n + 1) - \sigma, \\ \omega \tau_2 + \alpha &= \frac{\pi}{2} (4 n + 1) + \sigma. \end{aligned} \] Отметим, что $\tau_1 = (\tau_1^0, \tau_1^1, \dots)$ зависит от $n$.

    Аналогично находим точки переключения управления $u_2$: \[ \tag{10} \begin{aligned} \omega \tau_3 + \alpha &= \frac{\pi}{2} (4 n + 3) - \sigma, \\ \omega \tau_4 + \alpha &= \frac{\pi}{2} (4 n + 3) + \sigma, \end{aligned} \] где $\tau_3, \tau_4$ — точки включения и выключения соответственно.

    Так как вырожденный режим невозможен, то значение управления равно 0 или $\nu$.

    Заметим, что ширина ступеньки одна и та же (равна $2\sigma$). Надо найти $\alpha$ и $\sigma$. Начальная фаза находится из начальных условий, а вот $\sigma$ найти сложнее.

    Ищем $\alpha$. Рассмотрим общее условие трансверсальности: \[ \Delta \varphi_0 + \left.\paren{\lambda \Delta x - H \Delta t}\right|_{t_\text{н}}^{t_\text{к}} = 0. \] Из начальных условий следует, что \[ \varphi_0 = 0, \qquad \Delta x_\text{н} = \Delta x_\text{к} = \Delta t_\text{н} = 0. \] Тогда \[ H_\text{к} \Delta t_\text{к} = 0. \] Вариации незавимимы, поэтому $H_\text{к} = 0$. Заметим, что гамильтониан не зависит явно от $t$, то есть он является первым интегралом системы из уравнений движения и уравнений Эйлера-Лагранжа. Так как $H_\text{к} = 0$, а $H = \const$, то $H \equiv 0$.

    Рассмотрим $t_\text{н} = 0$: $u_1(0) = u_2(0) = 0$, поэтому \[ \begin{aligned} H &= \lambda_1 x_2(0) + \lambda (-\omega^2 x_1(0)) = \\ &= - a \omega x_2(0) \cos \alpha - a \omega^2 x_1(0) \sin\alpha = 0 \end{aligned} \] откуда \[ \tag{11} \tg \alpha = - \frac{x_2(0)}{\omega x_1(0)}. \]

    Осталость только найти $\sigma$ — займёмся этим позже.
  3. Фазовый портрет энергетически оптимальных траекторий.
    Введём обозначение: $u = u_1 - u_2$, причём в любой момент времени $u = \const$. Заметим, что оно может принимать только 3 значения: $(-\nu, 0, \nu)$. Тогда уравнение движения можно записать в виде \[ \ddot x_1 + \omega^2 x_1 = u. \] Решение этого уравнения выглядит в виде \[ \tag{12} \begin{aligned} x_1 &= b \sin(\omega t + \beta) + u \omega^{-2}, \\ x_2 &= b \omega \cos(\omega t + \beta), \end{aligned} \] где $b = \const \gt 0$ — амплитуда фазового движения, $\beta = \const$ — начальная фаза колебания. Какую траекторию определяет это движение?

    Так как $t_\text{н} = 0: u_1 = 0, u_2 = 0$, то \[ \tag{13} \tg \beta = \frac{\omega x_1(0)}{x_2(0)}, \qquad b = \sqrt{(x_1(0))^2 + (x_2(0))^2 \omega^{-2}} \] Отсюда следует, что \[ \begin{aligned} \sin \alpha &= - x_2(0), & \cos \alpha &= \omega x_1(0), \\ \sin \beta &= \omega x_1(0), & \cos \beta &= x_2(0). \end{aligned} \] То есть \[ \cos \alpha = \sin \beta, \qquad \sin \alpha = - \cos\beta, \] значит, фазы сдвинуты относительно друг друга на $\pi / 2$, то есть \[ \tag{14} \alpha = \beta \pm \frac{\pi}{2}. \] Так как частота одинаковая, то этот сдвиг по фазе сохраняется на всей траектории.

    Заметим, что $x_1$ сдвинут относительно $x_2$ на $\pi / 2$. Будем считать, что $\lambda_2$ и $x_2$ находится в противофазе.

    Для получения траектории исключим $t$. Получается, что траектории представляют собой окружности \[ \tag{15} \paren{\omega x_1 - u \omega^{-1}}^2 + x_2^2 = \omega^2 b^2 = \const. \] Оси: $\omega x_1, x_2$. Понятно, что в зависимости от управления будут получаться разные окружности.

    Построим линии переключения (тоже на рисунке). Так как \[ \begin{aligned} \omega \tau_1 + \alpha &= \frac{\pi}{2} - \sigma, \\ \omega \tau_2 + \alpha &= \frac{\pi}{2} + \sigma. \end{aligned} \] Отсюда найдём длительность активного участка: \[ \omega (\tau_2 - \tau_1) = 2 \sigma. \] Так как \[ \alpha = \beta - \frac{\pi}{2}, \] то \[ \begin{aligned} x_1 &= b \sin(\omega t + \beta) + u \omega^{-2} = \\ &= b \sin (\omega t + \alpha - \frac{\pi}{2}) + u \omega^{-2}. \end{aligned} \] Подставим сюда $t = \tau_1$, получаем \[ \begin{aligned} x_1 &= b \sin (\omega t + \alpha - \frac{\pi}{2}) + u \omega^{-2} = \\ &= b \sin ( \frac{\pi}{2} - \sigma + \frac{\pi}{2}) + u \omega^{-2} = \\ &= b \sin (\pi - \sigma) + u \omega^{-2} = \\ &= b \sin \sigma + u \omega^{-2}. \end{aligned} \] Аналогично \[ \begin{aligned} x_2 &= b \omega \cos (\omega \tau_1 + \beta) = \\ &= b \omega \cos (\pi - \sigma) = \\ &= -b \omega \cos \sigma. \end{aligned} \] Итак, \[ \begin{aligned} \omega x_1 &= \phantom{-} b \omega \sin \sigma + u \omega^{-1}, \\ x_2 &= - b \omega \cos \sigma. \end{aligned} \] Выразим $b \omega$ через $x_1$ и подставим в $x_2$: \[ x_2 = - (\omega x_1 - u \omega^{-1}) \ctg \sigma. \] Заметим, что $\ctg \sigma = \const$. Получили линейную зависимость, то есть прямую. Найдём наклон, взяв производную: \[ \dv{x_2}{(\omega x_1)} = - \ctg \sigma \bydef= \tg \delta, \] где $\delta$ — угол наклона. Отсюда следует, что \[ \delta = \frac{\pi}{2} + \sigma. \] Итак, на этой прямой лежат все точки включения управления $u_1$.

    Аналогично, подставив вместо $\tau_1$ точку $\tau_2$, получим, что угол наклона линии, на которой лежат точки выключения, равен \[ \delta = \frac{\pi}{2} - \sigma. \]

    Так как $\lambda_2$ находится с $x_2$ в противофазе, то (?).

    Сделать вывод относительно лучей, на которых эти точки переключения лежат, учитывая противофазу.

    Рассмотрим начальную точку $(\omega x_1^\text{н}, x_2^\text{н})$. Так как начинаем двигаться без управления, то точка движется по окружности с центром в точке $O$. После переключения движемся по окружности с центром в точке $O_1$. Промежуток активного движения длится $2 \sigma$. После выключения снова двигаемся по окружности с центром в точке $O$, после чего включается второе управление, поэтому начинаем двигаться по окружности с центром в $O_2$. И так далее.

    От чего зависит направление движения фазовой траектории? Вспомним, что $\dot x_1 = x_2$. В верхней полуплоскости (где $x_2 \gt 0$), точка будет двигаться слева направо ($\dot x_1 \gt 0$). В нижней полуплоскости соответственно $\dot x_1 \lt 0$, поэтому двигаемся справа налево.

    С каждым переключением управления меняется величина $b$, то есть радиус окружности. Определим, на сколько один промежуток управления уменьшает амплитуду колебания.

    2023-10-05

    Из геометрических соображений найдём, как меняется радиус после включения и выключения управления.

    Из треугольника $OO_1A$ по теореме косинусов \[ b_1^2 \omega^2 = b^2 \omega^2 + \nu^2 \omega^{-2} - 2 b \nu \cos\paren{\frac{\pi}{2} - \sigma}. \] Из треугольника $OO_1B$ по теореме косинусов \[ b_2^2 \omega^2 = b_1^2 \omega^2 + \nu^2 \omega^{-2} - 2 b_1 \nu \sin(2 \sigma + \gamma) \]

    Из треугольника $ACO_1$ получаем \[ \begin{aligned} \sin \gamma &= \frac{\nu \omega^{-1} - b \omega \sin \sigma}{b_1 \omega} \\ \cos \gamma &= \frac{b \omega \cos \sigma}{b_1 \omega}. \end{aligned} \] Тогда синус суммы углов равен \[ \begin{aligned} \sin (2 \sigma + \gamma) &= \frac{b}{b_1} \sin 2\sigma \cos \sigma - \frac{b}{b_1} \sin \sigma \cos 2\sigma + \frac{\nu}{b_1 \omega^2} \cos 2\sigma = \\ &= b_1^{-1} \paren{ b \sin\sigma + \nu \omega^{-2} \cos 2\sigma }. \end{aligned} \] Получим разницу между $b \omega$ и $b_1 \omega$: \[ \begin{aligned} b_2^2 \omega^2 &= b^2 \omega^2 + 2 \nu^2 \omega^{-2} - 4 b \nu \sin\sigma - 2 \nu^2 \omega^{-2} \cos 2\sigma = \\ &= (b \omega - 2 \nu \omega^{-1} \sin\sigma)^2, \end{aligned} \] откуда следует, что \[ b_2 = b - 2 \nu \omega^{-2} \sin \sigma, \tag{16} \] поэтому \[ b - b_2 = 2 \nu \omega^{-2} \sin \sigma. \tag{17} \] То есть за один "виток" управления амплитуда уменьшается на это число. Определим время, за которое мы сможем погасить колебание.

    За один период времени \[ T = \frac{2 \pi}{\omega} \] амплитуда уменшается на $4 \nu \omega^{-2} \sin \sigma$. Поэтому для погашения начальных колебаний амплитуды $b$ потребуется время \[ t_\text{д} = \frac{T b}{4 \nu \omega^{-2} \sin \sigma} \tag{18} \] (индекс "д" означает "действие").

    Рассчитаем функционал \[ J_1 = \int\limits_{t_\text{н}}^{t_\text{к}} (u_1 + u_2) dt. \] Управления включаются попеременно, причём каждая ступенька имеет ширину $2 \sigma$ и высоту $\nu$. Заметив, что включения и выключения происходят периодически, получаем, что \[ J_1 = 4 \sigma \nu \frac{t_\text{д}}{2 \pi} = \frac{4 \sigma \nu t_\text{д}}{\omega T}, \quad \text{где} \quad T = \frac{2 \pi}{\omega}. \] Откуда \[ J_1 = b \omega \frac{\sigma}{\sin \sigma} \tag{19} \] Так как \[ \lim_{\sigma \to 0} \frac{\sigma}{\sin \sigma} = 1, \] а при $\sigma \gt 0$ это отношение будет больше 1, то $J_1$ будет стремиться к минимуму, когда $\sigma \to 0$. Так как $J_1$ представляет собой сумму площадей под ступенями, то минимизация этих ступений как раз и является нашей задачей.

    Если $\sin \sigma \to 0$, то из формулы (18) $t_\text{д} \to \infty$. Это означает, что теоретически оптимальный режим такой, что кусочно-постоянное управление стремится к импульсному, а время действия стремится к бесконечности.

    Дальше считаем, что количество ступеней зафиксировано (значит, фиксировано и $t_\text{д}$). Теоретически можно считать, что количество положительных и отрицательных ступеней разное, хотя с точки зрения практики имеет смысл взять их одинаковыми.

  4. Нахождение ширины ступени управления

    Геометрически можно найти ширину ступеньки управления (см. Новосёлова, например), но мы построим общий алгоритм решения таких задач.

    Рассмотрим систему уравнений движения (6): \[ \begin{aligned} \dot x_1 &= x_2, \\ \dot x_2 &= - \omega^2 x_1 + u, \end{aligned} \] где \[ \tag{20} u = u_1 - u_2, \qquad x(0) = \paren{x_1^0, x_2^0}^T. \] В матричном виде эта система имеет вид \[ A = \begin{pmatrix} 0 & 1 \\ -\omega^2 & 0 \end{pmatrix}, \quad U = \paren{0, u}^T, \quad x_0 = \paren{x_1^0, x_2^0}^T = \paren{ \varphi_0, \dot\varphi_0 }^T. \]

    Если у $A$ вещественные собственные числа, то движение апериодическое.

    Матрица $A$ имеет пару чисто мнимых собственных значений $\pm i\omega$, следовательно, её жорданова форма $Y_A$ — диагональная матрица. Это значит, что существует такая неособая постоянная комплексная матрица $C$, которая преобразует $A$ к её жордановой форме: \[ \tag{21} C A C^{-1} = Y_A. \] Одна из матриц, удовлетворяющих (21), может быть такой: \[ C = \begin{pmatrix} \omega & -i \\ \omega & \phantom{-}i \end{pmatrix}. \] Сделаем замену переменных: $y = Cx$ тогда уравнения и условия (20) перейдут в следующие: \[ \tag{22} \dot y = Y_A y + V, \] где \[ y_0 = Cx_0, \qquad V = CU = \paren{ v_1, v_2 }^T = \paren{ -iu, iu }^T. \] Согласно структуре \[ Y_A = \begin{pmatrix} i \omega & 0 \\ 0 & -i \omega \end{pmatrix}, \] система (22) состоит из двух линейных неоднородных дифференциальных уравнений: \[ \tag{22*} \begin{aligned} \dot y_1 &= \phantom{-} i \omega y_1 + v_1, \\ \dot y_2 &= -i \omega y_2 + v_2. \end{aligned} \] Новые начальные условия $y_0 = Cx_0$ выглядят так: \[ \begin{aligned} y_1^0 &= \omega x_1^0 - i x_2^0, \\ y_2^0 &= \omega x_1^0 + i x_2^0. \end{aligned} \] Формализуем структуру функций управления через точки переключения: \[ \tag{23} u(t) = \nu \sum_{i=1}^{2r} (-1)^{i+1} H(t - \tau_i) + \nu \sum_{i=1}^{2q} (-1)^i H(t - \tilde \tau_i), \] где $p,q$ — количество положительных и отрицательный ступеней соответственно, а $H(t - \tau_i)$ — функция Хевисайда: \[ H(t - \tau_i) = \begin{cases} 1, & t - \tau_i \geqslant 0, \\ 0, & t - \tau_i \lt 0. \end{cases} \] Система (22) (или $(22^*)$) состоит из двух независимых друг от друга уравнений. Проинтегрируем их, учитывая при вычислении интегралов структуру управления $(23)$. Полученное решение системы (22) имеет вид \[ \begin{aligned} y_1 &= \left[ y_1^0 - \frac{\nu}{\omega} \paren{ F_1 - i F_2 } \right] \cos \omega t + \left[ \phantom{-}i y_1^0 - \frac{i \nu}{\omega} \paren{ F_1 - i F_2 } \right] \sin \omega t + \frac{u}{\omega}, \\ y_2 &= \left[ y_2^0 + \frac{\nu}{\omega} \paren{ F_1 - i F_2 } \right] \cos \omega t + \left[ -i y_2^0 - \frac{i \nu}{\omega} \paren{ F_1 + i F_2 } \right] \sin \omega t + \frac{u}{\omega}, \end{aligned} \] где \[ \begin{aligned} F_1 &= \sum\limits_{i=1}^{2r} (-1)^{i+1} H(t - \tau_i) \cos \omega \tau_i + \sum\limits_{i=1}^{2q} (-1)^i H(t - \tilde \tau_i) \cos \omega \tilde \tau_i, \\ F_2 &= \sum\limits_{i=1}^{2r} (-1)^{i+1} H(t - \tau_i) \sin \omega \tau_i + \sum\limits_{i=1}^{2q} (-1)^i H(t - \tilde \tau_i) \sin \omega \tilde \tau_i. \end{aligned} \] Чтобы решить поставленную задачу гашения колебания частоты $\omega$, потребуем, чтобы выражения в квадратных скобках (то есть при $\sin \omega t$ и $\cos \omega t$) обнулялись после отработки всех двигателей, то есть при $t \gt \tau_{2r}, t \gt \tau_{2q}$ \[ \tag{24} \begin{aligned} y_1^0 - \frac{\nu}{\omega} \paren{ F_1 - i F_2 } &= 0, & y_2^0 - \frac{\nu}{\omega} \paren{ F_1 + i F_2 } &= 0, \\ i y_1^0 - \frac{i \nu}{\omega} \paren{ F_1 - i F_2 } &= 0, & -i y_2^0 + \frac{i\nu}{\omega} \paren{ F_1 + i F_2 } &= 0. \end{aligned} \] Легко увидеть, что условия в каждом столбце выражений $(24)$ линейно зависимы друг от друга. Таким образом, из четырёх граничных условий остаётся только два, и в $(24)$ можно ограничиться только первой строкой.

    Формулы $(24)$ имеют комплексные коэффициенты, Чтобы перейти от них к вещественным коэффициентам, разделим каждое из условий $(24)$ на вещественные и мнимые части, принимая во внимание, что \[ y_1^0 = \omega x_1^0 - i x_2^0, \qquad y_2^0 = \omega x_1^0 + i x_2^0, \] тогда \[ \tag{24*} \left\{ \begin{aligned} \omega x_1^0 - \frac{\nu}{\omega} F_1 - i \paren{ x_2^0 - \frac{\nu}{\omega} F_2} &= 0, \\ \omega x_1^0 - \frac{\nu}{\omega} F_1 + i \paren{ x_2^0 - \frac{\nu}{\omega} F_2} &= 0. \end{aligned} \right. \] Чтобы условия $(24^*)$ выполнялись, необходимо и достаточно, чтобы их вещественная и мнимая части были равны нулю. Учитывая, что вещественная часть у двух равенств совпадает, а мнимые части отличаются лишь знаком, получим два новых граничных условия: \[ \tag{25} \omega x_1^0 - \frac{\nu}{\omega} F_1 = 0, \qquad x_2^0 - \frac{\nu}{\omega} F_2 = 0. \] Заметим, что все коэффициенты вещественны. С учётом $x_1^0 = \varphi_0, x_2^0 = \dot \varphi_0$, можно переписать условия $(25)$ в виде \[ K_1 = \omega \varphi_0 - \frac{\nu}{\omega} F_1 = 0, \qquad K_2 = \dot \varphi_0 - \frac{\nu}{\omega} F_2 = 0. \]

    Минимизируемый функционал \[ J_1 = \int\limits_{t_\text{н}}^{t_\text{к}} (u_1 + u_2) dt \] стоит переписать как функцию от точек переключения. Для этого учтём структуру управления, а также то, что и в начальный, и в конечный моменты времени управляющее воздействие не оказывается. Тогда \[ \tag{26} J_1 = \nu \left[ \sum\limits_{i=1}^{2r} (-1)^i \tau_i + \sum\limits_{i=1}^{2q} (-1)^i \tilde\tau_i \right]. \]

    расписать?
    Интегрируем по кускам, причём в начале каждого куска начальные условия — то состояние, какое получили на текущий момент.

    Далее будем решать задачу минимизации функционала $(26)$ при выполнении граничных условий $(25^*)$.

    Чтобы перейти к задаче безусловной минимизации, введём множители Лагранжа $\kappa_1$ и $\kappa_2$: \[ R = J_1 + \kappa_1 K_1 + \kappa_2 K_2. \] Согласно формулам $(25^*)$ и $(26)$, этот функционал зависит от точек переключения и множителей Лагранжа: \[ R = R(\kappa_1, \kappa_2, \tau_i, \tilde \tau_j), \qquad i = \overline{1, 2r}, \quad j = \overline{1, 2q}. \] Следовательно, необходимые условия минимума этой функции — условия экстремума функции многих переменных: \[ \tag{27} \pd{R}{\kappa_1} = 0, \quad \pd{R}{\kappa_2} = 0, \quad \pd{R}{\tau_i} = 0, \quad \pd{R}{\tilde\tau_j} = 0. \] Заметим, что согласно полученной нами структуре управления (см. рисунок 2) ступени управления чередуются и не пересекаются (введённая нами структура управления $u(t)$ явно эту информацию не содержит). К тому же, зная точки включения $\tau_1$ и выключения $\tau_2$ первого управления, все остальные точки переключения можно выразить через них. Эти два факта можно записать в виде \[ \tag{*} \begin{aligned} \tau_i &= \tilde \tau_i - \frac{\pi}{\omega} + \frac{2 \pi l}{\omega}, \\ \tau_i &= \tau_{i+2} - \frac{2 \pi}{\omega} + \frac{2 \pi l}{\omega}. \end{aligned} \] С учётом этого первые два условия $(27)$ могут быть переписаны в следующем виде: \[ \tag{28} \begin{aligned} K_1 &= \omega \varphi_0 + \frac{\nu}{\omega} (r + q) \left[ \cos \omega\tau_2 - \cos \omega\tau_1 \right] = 0, \\ K_2 &= \phantom{\omega} \varphi_0 + \frac{\nu}{\omega} (r + q) \left[ \sin \omega\tau_2 - \sin \omega\tau_1 \right] = 0. \end{aligned} \]

    Во втором уравнении стоит $\varphi_0$ или $\dot\varphi_0$?

    Заметим, что \[ \tau_2 - \tau_1 = 2 \sigma. \] Пусть $t^*$ — средний момент ступени управления, тогда \[ \tag{29} \tau_1 = t^* - \sigma, \qquad \tau_2 = t^* + \sigma. \] При подстановке формул $(29)$ в граничные условия $(28)$ приходим к следующим соотношениям: \[ \tag{30} \begin{aligned} K_1 &= \omega \varphi_0 - \frac{2}{\omega} \nu (r + q) \sin \omega t^* \sin \omega \sigma = 0, \\ K_1 &= \phantom{\omega} \dot\varphi_0 + \frac{2}{\omega} \nu (r + q) \cos \omega t^* \sin \omega \sigma = 0. \end{aligned} \] Функционал $J_1$ при подстановке $(*)$ и $(29)$ принимает вид \[ \tag{31} J_1 = 2 \nu (r + q) \sigma, \] что полностью соответствует полученной структуре управления (рис. 2) и тому смыслу, что функционал представляет собой суммарную площадь всех ступеней управления — математическая модель адекватна.

    Теперь $J_1, K_1, K_2$ — функции от $t^*$ и $\sigma$, а исходная задача оптимизации свелась к безусловной задаче минимизации функционала \[ R = J_1 + \kappa_1 K_1 + \kappa_2 K_2 \] по переменным $t^*, \sigma, \kappa_1, \kappa_2$: \[ \tag{27*} \pd{R}{t^*} = 0, \quad \pd{R}{\sigma} = 0, \quad \pd{R}{\kappa_1} = 0, \quad \pd{R}{\kappa_2} = 0. \] Два первых уравнения $(27^*)$ приводят к уравнениям \[ \tag{32} \begin{aligned} \paren{ \kappa_1 \sin \omega t^* + \kappa_2 \cos \omega t^* } \cos \omega \sigma &= 1, \\ \paren{ \kappa_1 \cos \omega t^* - \kappa_2 \sin \omega t^* } \sin \omega \sigma &= 0, \end{aligned} \] два последних же — к уравнениям $(30)$. Таким образом, уравнения $(30)$ и $(32)$ представляют собой систему из четырёх уравнений относительно четырёх неизвестных.

    Рассмотрим сначала уравнения $(32)$. Возможны два варианта:

    1. $\sin \omega\sigma = 0$. Тогда $\cos \omega\sigma = \pm 1$. В этом случае либо $\sigma = 0$ (управление отсутствует), либо $\sigma = \frac{\pi}{\omega}$, и управление непрерывно (см. $(*)$). Оба варианта не удовлетворяют требованиям поставленной задачи.
    2. $\sin \omega\sigma \neq 0$.

    Первая из формул $(32)$ даёт $\cos \omega\sigma \neq 0$, поэтому можно записать \[ \sin \omega\sigma = \sqrt{ 1 - \cos^2 \omega \sigma }. \]

    Согласно рис. 2 \[ \omega \sigma \lt \frac{\pi}{2}. \]
    Так как в этом случае \[ \kappa_1 \cos \omega t^* - \kappa_2 \sin \omega t^* = 0, \] то из первой формулы $(32)$ получаем \[ \begin{aligned} \sin \omega t^* &= \frac{\kappa_1}{\cos \omega \sigma (\kappa_1^2 + \kappa_2^2)}, \\ \cos \omega t^* &= \frac{\kappa_2}{\cos \omega \sigma (\kappa_1^2 + \kappa_2^2)}. \end{aligned} \] После подстановки этих выражений в $(30)$ получаем соотношение \[ \kappa_2 = - \frac{\dot \varphi_0}{\omega \varphi_0} \kappa_1. \] Учитывая его, из тех же граничных условий $(30)$ получаем уравнение для нахождения $\kappa_1$: \[ \paren{ \dot\varphi_0^2 + \omega^2 \varphi_0^2 } \kappa_1 = - \frac{2}{\omega} \nu (r + q) \sqrt{ \kappa_1^2 \paren{ \dot\varphi_0^2 + \omega^2 \varphi_0^2 } - \omega^2 \varphi_0^2 }. \] Отсюда \[ \tag{33} \kappa_1^2 = \frac {4 \nu^2 (r + q)^2 \omega^2 \varphi_0^2} {Z^2 \paren{ 4 \nu^2 (r + q)^2 - \omega^2 Z^2 }}, \quad \text{где} \quad Z^2 = \dot\varphi_0^2 + \omega^2 \varphi_0^2. \] Из второго уравнения $(32)$ следует, что если $\omega t^* \in [0; \pi/2) \cup (\pi/2; \pi]$, то $\omega t^*$ определяется однозначно: \[ \tg \omega t^* = \frac{\kappa_1}{\kappa_2} = - \frac{\omega \varphi_0}{\dot \varphi_0}. \] Тогда из первой формулы $(32)$ получаем выражение для $\omega \sigma$: \[ \tag{34} \cos \omega \sigma = \paren{ \kappa_1 \sin \omega t^* + \kappa_2 \cos \omega t^* }^{-1}, \] где:

    Так как $\omega \sigma \lt \dfrac{\pi}{2}$, то \[ \cos \omega \sigma \gt 0, \quad \sin \omega \sigma \gt 0. \] Следовательно, с учётом $(34)$ и $(35)$ выражения для половины ширины ступени управления окончательно выписываются в следующем виде: \[ \tag{36} \left\{ \begin{aligned} \cos \omega \sigma &= \frac{\omega z_0 \xi}{\abs{\kappa_1} Z}, \\ \sin \omega \sigma &= \frac {\sqrt{ Z^2 \kappa_1^2 - \omega^2 \varphi_0^2 }} {\abs{\kappa_1} Z}. \end{aligned} \right. \]

    Таким образом, ширину ступени управления $2 \sigma$ удалось получить с помощью явных формул $(36)$ через начальные условия и заданные параметры управления (максимальную мощность двигателей $\nu$ и количество включений управления $r+q$). Средний момент первой ступени задаётся формулами $(35)$

    Зная значения $t^*$ и $\sigma$, можно найти момент первого включения управления $\tau_1$ и момент выключения первой ступени управления $\tau_2$ по формулам $(29)$: \[ \tau_1 = t^* - \sigma, \qquad \tau_2 = t^* + \sigma. \] Остальные точки переключения $\tau_i$ и $\tilde \tau_i$ могут быть посчитаны с помощью $\tau_1$ и $\tau_2$ по формулам $(*)$: \[ \begin{aligned} \tau_i &= \tilde \tau_i - \frac{\pi}{\omega} + \frac{2 \pi l}{\omega}, \\ \tau_i &= \tau_{i+2} - \frac{2 \pi}{\omega} + \frac{2 \pi l}{\omega}, & l = 0, 1, \dots \end{aligned} \]

§2. Оптимальное гашение плоских колебаний спутника по быстродействию

Пусть $u = u_1 - u_2$, тогда $-\nu \leqslant u \leqslant \nu$. Рассматриваем задачу быстродействия: \[ \int\limits_{t_\text{н}}^{t_\text{к}} dt, \implies F_0 = 1. \] Составим гамильтониан: \[ H = - F_0 + \lambda f = -1 + \lambda f \sim \lambda f. \] Так как \[ \begin{aligned} \dot x_1 &= x_2, \\ \dot x_2 &= - \omega^2 x_1 + u, \end{aligned} \] то \[ H = \lambda_1 x_2 + \lambda_2 (- \omega^2 x_1 + u). \] Гамильтониан линейно зависит от управление, следовательно, из ПМП следует, что \[ u = \nu \cdot \sign \lambda_2(t). \] Составим уравнения Эйлера-Лагранжа: \[ \dot \lambda_i = - \pd{H}{x_i}, \] то есть \[ \begin{aligned} \dot \lambda_1 &= \omega^2 \lambda_2, \\ \dot \lambda_2 &= - \lambda_1, \end{aligned} \] откуда можно найти решение \[ \begin{aligned} \lambda_1 &= -a \omega \cos(\omega t + \alpha), \\ \lambda_2 &= a \sin(\omega t + \alpha), \end{aligned} \] где $a = \const \gt 0$ и $\alpha = \const$.

Вырожденный режим \[ \pd{H}{u_1} = \pd{H}{u_2} = 0 \] невозможен, поэтому \[ \pd{H}{u} = \lambda_2 = 0 \] задаёт точки переключения.

Для получения траектории решим исходное уравнение и исключим время: \[ \paren{\omega x_1 - u \omega^{-1}}^2 + x_2^2 = \omega^2 b^2 = \const. \] Фазовый портрет:

Пусть $x_0 = (1, 1), t_0 = 0, \omega = 1, \nu = 1$. Тогда \[ O_1(1,0), \qquad O_2(1, 0). \]
Из геометрических соображений: \[ \tg \varphi = \frac{1}{2}, \qquad \sin \varphi = \frac{1}{\sqrt{5}} \qquad \cos \varphi = \frac{2}{\sqrt{5}}. \] Значит, \[ \sin 2 \varphi = \frac{4}{5}, \implies 2 \varphi = \arcsin \frac{4}{5} \approx 0.3 \pi. \] Так как $\varphi = \omega t$, а $\omega = 1$, то $t = \varphi$. Значит, \[ t_\text{к} = 2 \varphi + \frac{\pi}{2} \approx 0.3\pi + 0.5 \pi = 0.8\pi. \] Из уравнений Эйлера-Лагранжа известно, что \[ \lambda_2 = a \sin(t + \alpha). \] Найдём точки переключения управления. Для этого решим уравнение $\lambda_2 = 0$: \[ a \sin(\tau + \alpha) = 0. \] Так как из рисунка следует, что $\tau = 2 \varphi$, то $\alpha = - 2 \varphi$. Значит, \[ u = \nu \cdot \sign \lambda_2 \implies u = \sign \sin(t - 2 \varphi). \]

2023-10-?

Я запутался в датах... В общем, по темам эта лекция стоит тут.

Глава 3. Управление живыми системами

§1. Математические модели экологии

1. Уравнение отдельной популяции
  1. Пусть $N(t)$ — численность населения, а $\lambda = \const$ — коэффициент, означающий разницу между рождаемостью и смертностью. Тогда \[ \dot N(t) = \lambda N(t), \qquad N(0) = N_0 \gt 0 \] называют моделью Мальтуса. Минус этой модели заключается в том, что она не учитывает ограниченность ресурсов.
  2. Пусть $K$ — ёмкость среды, тогда \[ \dot N(t) = \lambda \paren{ 1 - \frac{N(t)}{K} } N(t), \qquad t \geqslant 0, \quad N(0) = N_0 \gt 0 \] называют уравнением Ферхюльста (или логистическим уравнением). Коэффициент $\dfrac{\lambda}{K}$ называют коэффициентом внутривидовой борьбы.

Рассмотрим подробнее логистическое уравнение. Оно является уравнением с разделяющимися переменными: \[ \frac{K dN}{N K - N^2} = \lambda dt, \] откуда можно найти, что \[ \tag{3} N(t) = \lambda N_0 e^{\lambda t} {\left[ \lambda + \frac{\lambda}{K} N_0 (e^{\lambda t} - 1) \right]}^{-1} = \frac{\lambda N_0 K e^{\lambda t}} {\lambda K + \lambda N_0 e^{\lambda t} - \lambda N_0}. \] Это уравнение называют логистической кривой.

Покажем, что решение логистического уравнения асимптотически стремится к ёмкости среды. Действительно, \[ \begin{aligned} \lim\limits_{t \to \infty} \frac{1}{N(t)} &= \lim\limits_{t \to \infty} \frac{\lambda K + \lambda N_0 e^{\lambda t} - \lambda N_0} {\lambda N_0 K e^{\lambda t}} = \\ &= \lim\limits_{t \to \infty} \left[ \frac{1}{N_0 e^{\lambda t}} + \frac{1}{K} - \frac{1}{K e^{\lambda t}} \right] = \\ &= \frac{1}{K}, \end{aligned} \] откуда и следует, что \[ \lim\limits_{t \to \infty} N(t) = K. \]

Найдём теперь положение равновесия логистического уравнения. Для этого приравняем нулю его правую часть: \[ \lambda \paren{ 1 - \frac{N(t)}{K} } N(t) = 0, \] откуда следует, что \[ N(t) \equiv N_1 = 0 \quad \text{и} \quad N(t) \equiv N_2 = K \] являются положениями равновесия.

Рассмотрим систему в отклонениях. Для этого введём вспомогательную переменную \[ \xi = N(t) - N_i \implies N(t) = \xi + N_i, \] тогда \[ \begin{aligned} \dot \xi &= \lambda \paren{ 1 - \frac{\xi + N_i}{K} } (\xi + N_i) = \\ &= \lambda (\xi + N_i) - \frac{\lambda}{K} \paren{\xi + N_i}^2 = \\ &= \lambda (\xi + N_i) - \frac{\lambda}{K} \paren{\xi^2 + 2 \xi N_i + N_i^2}. \end{aligned} \] Отбрасывая члены порядка малости выше 1, окончательно получаем линеаризованное уравнение \[ \dot \xi = \underbrace{\paren{\lambda - 2 N_i \frac{\lambda}{K}}}_{f'(N_i)} \xi + \underbrace{\lambda N_i - \frac{\lambda}{K} N_i^2}_{\const}. \] Итак, получаем два случая:

  1. Если $N_i = N_1 = 0$, то $f'(N_i) = \lambda \gt 0$, то есть положение $N(t) \equiv 0$ неустойчиво.
  2. Если $N_i = N_2 = K$, то $f'(N_i) = - \lambda \lt 0$, то есть положение $N(t) \equiv K$ асимптотически устойчиво. Отсюда следует, что если $N_0 \gt K$, то $\dot N \lt 0$, а также если $N_0 \lt K$, то $\dot N \gt 0$.

Обобщение логистического уравнения — модель с запаздыванием Хатчинсона: \[ \dot N(t) = \lambda \paren{ 1 - \frac{N(t - h)}{K} } N(t), \qquad h \gt 0. \]

Пусть $N(t, \tau)$ — та часть популяции, возраст особей в которой не превышает $\tau$. Тогда \[ \pd{N(t, \tau)}{\tau} = n(t, \tau), \] где $n(t, \tau)$ — плотность. Тогда \[ \pd{n(t, \tau)}{t} + \pd{n(t, \tau)}{\tau} = - \mu(t, \tau) n(t, \tau) + g(t, \tau), \] где $\mu(t, \tau) n(t, \tau)$ отвечает за смертность, а $g(t, \tau)$ — за миграцию.

Начальные условия: \[ n(0, \tau) = n_0(\tau), \qquad \tau \geqslant 0. \]

Можем управлять этой системой, влияя на рождаемость: \[ n(t, 0) = u(t) \int\limits_{\tau_1}^{\tau_2} k(t, \tau) n(t, \tau) d\tau, \] где $\tau_1, \tau_2$ — границы репродуктивного возраста, а $k(t, \tau)$ — вклад в рождаемость.

Если учитывать неоднородность по территории, то можно ввести в рассмотрения обознчаения: Модель: \[ \pd{D(t, x, y)}{t} = \alpha \left[ \pdv2{D(t, x, y)}{x} + \pdv2{D(t, x, y)}{y} \right] + \lambda D(t, x, y) \paren{ 1 - \frac{D(t - h, x, y)}{K} }. \]
2. Сообщества двух и более видов

Классическая модель — модель Лотки-Вольтерры («хищник-жертва»): \[ \begin{aligned} \dot x(t) &= a_1 x - a_2 x y, \\ \dot y(t) &= a_3 a_2 x y - a_4 y, \end{aligned} \] где

Приравнивая правые части нулю, получаем стационарное положение: \[ \tilde x = \frac{a_4}{a_2 a_3}, \qquad \tilde y = \frac{a_1}{a_2}. \]

Домножив первое уравнение на $(-a_2 a_3)$, второе — на $(-a_2)$ и сложив полученные уравнения, получим \[ a_2 a_3 \dot x + a_2 \dot y = a_2 a_3 a_1 x - a_2 a_4 y. \] С другой стороны, разделив первое на $(a_4 x)$, второе на $(a_1 y)$ и сложив, получим \[ \frac{a_4 \dot x}{x} + \frac{a_1 \dot y}{y} = - a_2 a_4 y + a_2 a_3 a_1 x. \]

Заметим, что у полученных уравнений совпадают правые части, поэтому \[ \frac{a_4 \dot x}{x} + \frac{a_1 \dot y}{y} = a_2 a_3 \dot x + a_2 \dot y. \] Решая это дифференциальное уравнение, приходим к неявной форме записи решения: \[ a_2 a_3 x + a_2 y = a_4 \ln x + a_1 \ln y + \ln C, \] или \[ e^{a_2 a_3 x} x^{-a_4} = C e^{-a_2 y} y^{a_1}. \]

Если сделать замену \[ x_1(t) = \frac{x(t)}{\tilde x}, \qquad y_1(t) = \frac{y(t)}{\tilde y}, \] то $\tilde x_1 = \tilde y_1 = 1$, и уравнение траектории упростится: \[ \paren{ \frac{e^{x_1}}{x_1} }^{a_4} \paren{ \frac{e^{y_1}}{y_1} }^{a_1} = C_1. \]

Выясним теперь период. Рассмотрим малые отклонения: \[ \xi(t) = x(t) - \tilde x, \qquad \eta(t) = y(t) - \tilde y, \] тогда \[ \begin{aligned} \dot \xi(t) &= -\frac{a_4}{a_3} \eta(t), \\ \dot \eta(t) &= a_1 a_3 \xi(t). \end{aligned} \] Определитель этой системы: \[ \begin{vmatrix} -\lambda & - \frac{a_4}{a_3} \\ a_1 a_3 & - \lambda \end{vmatrix}, \] откуда $\lambda_{1,2} = \pm i \sqrt{a_1 a_4}$. Вводя обозначение $\omega = \sqrt{a_1 a_4}$, можно будет записать период как \[ T = \frac{2 \pi}{\sqrt{a_1 a_4}} = \frac{2 \pi}{\omega}. \]

Можно учесть внутривидовую конкуренцию, тогда модель Лотки-Вольтерры преобразуется в следующую модель: \[ \begin{aligned} \dot x(t) &= a_1 x - a_2 x y - \gamma_1 x^2, \\ \dot y(t) &= a_3 a_2 x y - a_4 y - \gamma_2 y^2. \end{aligned} \]
Пусть: Тогда модель можно записать так: \[ \begin{aligned} \dot N_i(t) &= N_i(t) \left[ \lambda_i - \sum\limits_{j=1}^{n} \alpha_{ij} N_j(t) - \sum\limits_{j=1}^{m} \beta_{ij} P_j(t) \right], & i &= \overline{1,n}, \\ \dot P_p(t) &= P_p(t) \left[ -\gamma_p + \sum\limits_{j=1}^{m} \delta_{pj} N_j(t) \right], & p &= \overline{1,p}, \end{aligned} \] где $\lambda_i, \alpha_{ij}, \beta_{ij}, \gamma_p, \delta_{pj} = \const \gt 0$.

2023-10-26

§2. Задача управления в экологии

1. Постановка задачи управления экосистемами
Важные замечания:
  1. Управление однонаправленное (внесли удобрение — уже не выкопаем).
  2. Фазовые переменные всегда больше нуля, иначе теряется смысл.
  3. Управление может быть либо прямым (убил зайца), либо опосредованным (убил зайца — повлиял на популяцию волков).
Примеры задач:
  1. Оптимальный вылов популяции.

    Пусть $u(t)$ — интенсивность вылова (количество выловленных особей за единицу времени). Промежуток времени $t \in [0; T]$. Требуется рассчитать максимальный вылов.

    Уравнение Ферхюльста: \[ \tag{1} \dot N(t) = \lambda N(t) \paren{ 1 - \frac{N(t)}{K} } - u (t) N(t), \] где $K$ — ёмкость среды.

    Критерий качества: \[ \int\limits_{0}^{T} u(t) N(t) dt \to \max \]

    Ограничения: \[ 0 \leqslant u \leqslant \gamma = \const. \]

    Если $N_0$ — желаемое количество особей в оставшейся популяции, то критерий качества усложняется: \[ \int\limits_{0}^{T} \left[ \alpha_1 u(t) N(t) - \alpha_2 \paren{N(t) - N_0}^2 \right] dt \to \max, \] где $\alpha_1, \alpha_2 \gt 0$ — некоторые константы (весовые коэффициенты), которые влияют на предпочтительность того или иного критерия (условно: количество выловленной рыбы или остаток в популяции).

  2. Управление численностью народонаселения.

    С прошлого занятия:

    Пусть $N(t, \tau)$ — численность, где $\tau$ — крайний возраст в популяции. Пусть $n(t, \tau)$ — плотность: \[ \pd{N(t, \tau)}{\tau} = n(t, \tau). \] Тогда можно рассмотреть модель \[ \tag{5} \pd{n(t, \tau)}{t} + \pd{n(t, \tau)}{\tau} = - \mu(t, \tau) n(t, \tau) + g(t, \tau), \] где $\mu(t, \tau)$ — смертность, $g(t, \tau)$ — миграция.

    Управляющий параметр — рождаемость.

    \[ \tag{6} n(t, 0) = u(t) \int\limits_{\tau_1}^{\tau_2} K(t, \tau) n(t, \tau) d\tau, \] где $u(t)$ — коэффициент скорости рождаемости, $K(t, \tau)$ — вклад каждой возрастной группы, $n(t, 0)$ — родившиеся в момент времени $t$. \[ \tag{7} n(0, \tau) = n_0(\tau), \qquad \tau \geqslant 0. \]
    Определим критерий качества: \[ \int\limits_{0}^{\infty} d\tau \int\limits_{0}^{T} \left[ n(t, \tau) - n_0(\tau) \right]^2 h(t, \tau) dt \to \inf\limits_u, \] где: Ограничения: \[ 0 \leqslant u(t) \leqslant u_0 = \const. \] Управление $u(t)$ ходит только в граничное условие (видимо, $(6)$). В $(5)$ — возможный вариант управления: $u(t) = g(t, \tau)$.
  3. Управления системой хищник-жертва.

    Рассматриваем уравнения Лотки-Волтерры. Эффективность управления жертвами: \[ r_1 u(t) x(t), \] эффективность управления хищниками: \[ r_2 u(t) y(t). \] Например, $u(t)$ может быть интенсивностью уничтожения, например, а $r_1, r_2 \geqslant 0$ — заданные константы.

    \[ \tag{1} \begin{aligned} \dot x(t) &= x(t) \left[ a_1 - a_2 y(t) - r_1 u(t) \right], \\ \dot y(t) &= y(t) \left[ a_3 a_2 x(t) - a_4 - r_2 u(t) \right], \end{aligned} \] с начальными условиями \[ \tag{2} x(0) = x_0 \gt 0, \qquad y(0) = y_0 \gt 0. \] Ограничение на управление: \[ \tag{3} 0 \leqslant u(t) \leqslant \gamma = \const. \] Если мы находимся в положении экологического равновесия $R(\tilde x, \tilde y)$, то цели управления могут быть разные. Основных две:

Управление по быстродействию системой хищник-жертва

Будем рассматривать задачу (1)-(3). Положим $r_2 = 0$, то есть расматриваем опосредованное управление хищниками через популяцию жертв.

Приравняв правые части к нулю, получим положение равновесия (ненулевое). \[ \tag{4} \tilde x = \frac{a_4}{a_2 a_3}, \qquad \tilde y = \frac{a_1}{a_2}. \] Цель: за наикратчайшее время привести систему (1) из положения (2) в положение (4).

Критерий качества: \[ \tag{5} J = \int\limits_{0}^{T} dt \to \min. \] Граничные условия: \[ x(T) = \tilde x, \qquad y(T) = \tilde y. \]

Введём новые переменные: \[ \begin{aligned} t &= \frac{t'}{a_1}, & x &= \frac{a_4 x'}{a_2 a_3}, & y &= \frac{a_1 y'}{a_2}. \\ b &= \frac{a_4}{a_1}, & \gamma &= \frac{\gamma' r_1}{a_1}, & u &= \frac{a_1 u'}{r_1}. \end{aligned} \] Подставляем всё в уравнение (1). Чтобы не нагромождать, уберём штрихи.

Получаем систему в новых переменных: \[ \tag{6} \begin{aligned} \dot x(t) &= (1 - y - u) x, \\ \dot y(t) &= b y (x - 1). \end{aligned} \] Начальные условия: \[ x(0) = x_0 \gt 0, \qquad y(0) = y_0 \gt 0, \] ограничения: \[ 0 \leqslant u \leqslant \gamma, \qquad \gamma \gt 0, \quad t \gt 0. \] Тогда положение равновесия представляется как $R(1, 1)$.

Задача: перевести систему из $(x_0, y_0)$ в положение $R(1, 1)$ за кратчайшее время.

Заметим, что задача стационарна (правые части от $t$ не зависят) с незаданным временем. В этом случае функция Гамильтона составляется как \[ H(t, x, y, \lambda, u) = \lambda(t) f + \lambda_0 F_0, \] причём \[ \begin{gathered} \lambda_0 \leqslant 0, \qquad \lambda_0 \equiv \const, \\ \max H = 0. \end{gathered} \]

Тогда \[ \tag{7} H = \lambda_1 \paren{ 1 - y - u } x + \lambda_2 b y (x - 1) + \lambda_0, \] где $\lambda_1, \lambda_2$ удовлетворяют уравнениям Эйлера-Лагранжа: \[ \tag{8} \begin{aligned} \dot \lambda_1 &= - \pd{H}{x} = - \lambda_1(t) (1 - y - u) - \lambda_2(t) b y, \\ \dot \lambda_2 &= - \pd{H}{y} = \lambda_1(t) x + \lambda_2(t) b (1 - x). \end{aligned} \] Тогда \[ \tag{9} \max\limits_{0 \leqslant u \leqslant \gamma} H(x, y, \lambda_0, \lambda_1, \lambda_2, u) = H(x, y, \lambda_0, \lambda_1, \lambda_2, u_0) = 0. \] Перегруппируем функцию Гамильтона: \[ H = \lambda_1 (1 - y) x + \lambda_2 b y (x - 1) + \lambda_0 - \lambda_1 u x. \] Отсюда следует, что \[ \max\limits_u H = \max\limits_u \paren{ - \lambda_1 u x }. \] Учитывая, что $x \gt 0$, исключаем его из исследования. Получаем, что всё зависит от того, какой вид имеет $\lambda_1$, поэтому \[ u = \begin{cases} 0, & \lambda_1 \geqslant 0, \\ \gamma, & \lambda_1 \lt 0. \end{cases} \] Определим, сколько раз и в какие моменты $\lambda_1(t)$ меняет знак, откуда мы получим условия на изменения $u(t)$ во времени.

Для исследования того, когда же $\lambda_1(t)$ меняет знак, перейдём к новым переменным: \[ \tag{*} \varphi_1(t) = x(t) \lambda_1(t), \qquad \varphi_2(t) = y(t) \lambda_2(t). \] Заметим, что знаки $\varphi_1, \varphi_2$ совпадают со знаками $\lambda_1, \lambda_2$, поэтому будем исследовать поведение $\varphi_1(t)$.

Возьмём производную по времени от $\varphi_1(t),\varphi_2(t)$ в силу системы \[ \tag{10} \begin{aligned} \dot \varphi_1(t) &= - b x(t) \varphi_2(t), \\ \dot \varphi_2(t) &= \phantom{-} y(t) \varphi_1(t). \end{aligned} \] Из (6) видим, что знак $\dot y(t)$ зависит от того, как ведёт себя $x$, поэтому надо рассмотреть два случая: $x \gt 1$ и $x \lt 1$.

  1. Пусть $x(t) \gt 1$ при $t \in [t_1, t_2]$, тогда $\dot y(t) \gt 0$. Тогда в момент переключения управления $\tau \in [t_1, t_2]$ \[ \pd{H}{u} = 0, \implies \pd{H}{u} = \lambda_1(\tau) x(\tau) = \varphi_1(\tau) = 0. \] Из уравнений (6) и (9) следует, что \[ \tag{11} \frac{\varphi_1(t) \dot x(t)}{x(t)} + \frac{\varphi_2(t) \dot y(t)}{y(t)} + \lambda_0 = 0. \] Тогда \[ \varphi_2(t) = \frac{y}{\dot y} \paren{ - \lambda_0 - \varphi_1(t) \frac{\dot x}{x} }, \qquad t \in [t_1, t_2]. \] Подставим это выражение в (10), получим дифференциальное уравнение относительно $\varphi_1(t)$: \[ \tag{12} \begin{aligned} \dot \varphi_1(t) &= - b x(t) \frac{y}{\dot y} \paren{ - \lambda_0 - \varphi_1(t) \frac{\dot x}{x} }, \\ \varphi_1(\tau) &= 0. \end{aligned} \] Получили задачу Коши. Уравнение неоднородное. Известно, что решение уравнения \[ \dot z(t) = p(t) z(t) + f(t) \] представляется в виде \[ z(t) = e^{-\int\limits_{t_0}^{t} p(t) dt} \paren{ z_0 + \int\limits_{t_0}^{t} f(s) e^{\int\limits_{t_0}^{t} p(s) ds} ds }. \] Тогда знак функции определяется знаком $f(s)$. В нашем случае неоднородность: \[ f(t) = b x(t) \lambda_0 \frac{y}{\dot y}. \] Известно, что $b, x, \dot y, y \gt 0$, а $\lambda_0 \leqslant 0$.

    Покажем, что $\lambda_0 \neq 0$. От противного: пусть $\lambda_0 = 0$, тогда $f(t) \equiv 0$, поэтому \[ \dot \varphi_1(t) = \frac{b}{\dot y} y \dot x \varphi_1. \] Решение: \[ \varphi_1(t) = \varphi_1(t_0) e^{- \int\limits_{t_0}^{t} b y \dot x (\dot y)^{-1} ds}. \] В нашем случае $t_0 = \tau$, поэтому \[ \varphi_1(t) \equiv 0. \] Раз это так, то из (*) следует, что $\lambda_1(t) \equiv 0$, а поэтому и $\dot \lambda_1(t) \equiv 0$. Но тогда из (8) следует, что $\lambda_2(t) \equiv 0$. Но по условию принципа максимума $\lambda(t) \not \equiv 0$ — пришли к противоречию. Отсюда следует, что $\lambda_0 \lt 0$.

    Значит, $f(\tau) \lt 0$, откуда следует, что $\varphi_1(t) \lt 0$ при $t \in [\tau, t_2]$. Тогда в обратном интегрировании (от $\tau$ до $t_1$) $\varphi_1(t) \gt 0$.

    Итак, \[ \varphi_1(t) = \begin{cases} \gt 0, & t \in [t_1, \tau), \\ = 0, & t = \tau, \\ \lt 0, & t \in (\tau, t_2]. \end{cases} \] Отсюда следует, что \[ u_0 = \begin{cases} 0 & t \in [t_1, \tau], \\ \gamma, & t \in (\tau, t_2]. \end{cases} \]

  2. Пусть теперь $x(t) \lt 1$. Аналогичным образом выясняем, что \[ u_0 = \begin{cases} \gamma & t \in [t_1, \tau], \\ 0, & t \in (\tau, t_2]. \end{cases} \]

Как выглядит фазовая траектория? Надо найти линии переключения. Их строим двумя частями: аналитически и численно. Выясняется, что вид линий переключения зависит от $\gamma$.

  1. $\gamma = 1$
    $APRSB$ — линии переключения, которые строятся двумя способами: $APR$ — аналитически, а $RSB$ — численно.

    Рассмотрим второе уравнение системы (6): \[ \dot y = b y(x - 1). \] Если $x \gt 1$, то $\dot y \gt 0$, поэтому движение будет восходящим справа от точки равновесия. Если же $x \lt 1$, то $\dot y \lt 0$, поэтому движение будет нисходящим. Получается, что движение происходит против часовой стрелки.

    Алгоритм получения $RSB$ одинаков для всех $\gamma$. Согласно линии переключения в случае, если $x \gt 1$, $u = \gamma = 1$, и $u = 0$ иначе (см. рисунок).

    Для численного построения выходим из $R(1, 1)$ и начинаем интегрировать при $t = -t$. Тогда уравнения движения имеют вид \[ \tag{13} \begin{aligned} \dot x &= (u + y - 1) x(t), \\ \dot y &= b (1 - x) y(t). \end{aligned} \] Также \[ \tag{14} \begin{aligned} \dot \varphi_1 &= b x \varphi_2, \\ \dot \varphi_2 &= - y \varphi_1. \end{aligned} \] Так как справа от $R$ управление $u = \gamma$, то $\varphi_1(0) \lt 0$. Тогда \[ \tag{15} x(0) = 1, \qquad y(0) = 1, \qquad \varphi_1(0) \lt 0, \qquad \varphi_2(0) \text{ — любое}. \] Интегрируем (13) и (14) при начальных условиях (15). Доходим то точки, когда $\varphi_1(\tau_1) = 0$, переключаем управление на $u = 0$. Считаем \[ \tag{15'} x(\tau_1), \qquad y(\tau_1), \qquad \varphi_2(\tau_1) \] и решаем уравнения (13)-(14) при (15'). Доходим до точки, где $\varphi_1(\tau_2) = 0$. Эта точка как раз и принадлежит $RSB$.

    Придавая различные значения $\varphi_1(0), \varphi_2(0)$, получаем $RSB$.


    $APR$ является траекторией решения уравнений (13) при $u = 0$: \[ \tag{16} \begin{aligned} \dot x &= xy, \\ \dot y &= b (1 - x) y \end{aligned} \] с начальными условиями \[ x(0) = 1, \qquad y(0) = 1. \] Интеграл имеет вид \[ \paren{x e^{-x}}^b e^{-y} = C. \] Найдём $C$, при начальных условиях \[ C = e^{-b-1}, \] тогда уравнение траектории (линии $APR$): \[ \ln x^b = xb + y - b - 1. \] Можно выразить $y$: \[ \tag{17} y = b \ln x - xb + b + 1. \] Из (16) $x \gt 0$, поэтому $\dot x \gt 0$. Так как $x \gt 1$, то $\dot y \lt 0$. Если предположим, что $x \limto{t \to \infty} \infty$, то из (17) получится, что $y \limto{t \to \infty} \lt 0$. Отсюда следует, что \[ x \limto{t \to \infty} x^* \lt \infty. \] Если бы \[ y \limto{t \to \infty} a, \] то из (16) получилось бы, что $\dot x = ax$, то есть $x = e^{ax}$, откуда следовало бы, что \[ x \limto{t \to \infty} \infty, \] что невозможно. Отсюда следует, что $y \limto{t \to \infty} 0$.

  2. $\gamma \gt 1$
    Аналогично
  3. $\gamma \lt 1$ — более сложная ситуация.
    Две части строятся численно: $NA$ и $RSB$. Точек переключения становится больше, спираль накручивается. Это следует из того, что интенсивность управления низкая, поэтому точек переключения много.

2023-11-02

Модели культивирования микроорганизмов и управление ими

Непрерывная культура микроорганизмов. Модель Мано

Самый распространённый культиватор — хемостат. Пусть объём культиватора $V$, скорость притока субстрата $f$, тогда скорость протока \[ D = \frac{f}{V}. \] Если она больше максимальной, то культура вымывается — естественное ограничение.

Пусть $x$ — плотность биомассы (культуры), $s$ — концентрация субстрата. Считаем, что всё, что находится в культиваторе, равномерно (плотность и культуры, и субстрата в каждой точке одинаковы), а процесс непрерывный.

Модель Мано: \[ \tag{1} \begin{aligned} \dot x &= \mu(s) x - Dx, \\ \dot s &= -\alpha \mu(s) x + D (s_0 - s), \\ \mu(s) &= \frac{\mu_m s}{K_s + s}, \end{aligned} \] где:

Для упрощения модели введём переменные \[ x' = \frac{\alpha x}{K_s}, \qquad y = \frac{s}{K_s}, \qquad y_0 = \frac{s_0}{K_s}, \qquad t' = t \mu_m, \qquad D' = \frac{D}{\mu_m}. \] Тогда (опускаем штрихи) \[ \tag{2} \begin{aligned} \dot x &= \phantom{-} \mu(y) x - Dx, \\ \dot y &= - \mu(y) x + D (y_0 - y), \\ \mu(y) &= \frac{y}{1 + y}. \end{aligned} \]

Найдём стационарные положения (приравниваем правые части к нулю). Получаем два положения: \[ \tag{3} \tilde x_1 = 0, \qquad \tilde y_1 = y_0 \] и \[ \tag{4} \tilde x_2 = y_0 - \frac{D}{1 - D}, \qquad \tilde y_2 = \frac{D}{1 - D}. \]

Положение (4) — единственное положение, которое подходит по практическим соображениям (в (3) нет культуры).

В силу ограничения \[ y_0 = \frac{s_0}{K_s} \] положение (4) имеет смысл только при \[ D \leqslant \frac{y_0}{1 + y_0} = D_B, \] где $D_B$ — скорость вымывания. Заметим, что $D_B \lt 1$. В исходных переменных \[ D_B = \frac{\mu_m s_0}{K_s + s_0}. \]

Введём переменные \[ \tag{6} \xi = x - \tilde x_2, \qquad \eta = y - \tilde y_2. \] Подставим их в уравнение (2), а полученные уравнения (в малых отклонениях) линеаризуем около положения равновесия. Тогда определитель линеаризованной системы будет иметь вид \[ \begin{vmatrix} - \lambda & (D_B - D) (1 + y_0) (1 - D) \\ -D & -(D_B - D) (1 + y_0) (1 - D) - D - \lambda \end{vmatrix} = 0. \] Решение этого характеристического уравнения даёт два решения: \[ \begin{aligned} \lambda_1 &= -D, \\ \lambda_2 &= -(D_B - D)(1 + y_0) (1 - D). \end{aligned} \] Получаем, что (4) — асимптотически устойчивое положение равновесия (т.е. рабочее состояние культиватора).


Задачу максимизации выхода биомассы можно найти в тимсе. Она решается геометрическими методами (ПМП не нужен).

Минимизация времени выхода культиватора на стационарный режим

Будем рассматривать модель (2), но для наглядности вместо $y$ будем использовать $s$: \[ \tag{7} \begin{aligned} \dot x &= \mu(s) x - D x, \\ \dot s &= - \mu(s) x - Ds + D s_0, \\ \mu(s) &= \frac{s}{1 + s}. \end{aligned} \] В качестве управления будем рассматривать концентрацию подаваемого субстрата $s_0$. Ограничения: \[ 0 \leqslant s_0 \leqslant a. \] Будем рассматривать стационарное положение \[ \tag{8} \tilde x = s_0 - \frac{D}{1 - D}, \qquad \tilde s = \frac{D}{1 - D}. \] Надо построить $s_0 = s_0(x, s)$ переводящее систему (7) в состояние (8) за наименьшее время.

Обозначим через $s_p$ выходную концентрацию субстрата. Ясно, что $s_p \lt a$. Управление представимо в виде \[ s_0 = s_p + u, \] где \[ -s_p \leqslant u \leqslant (a - s_p) =: u_0, \qquad u_0 = \const, \; u_0 \gt 0. \] Будем теперь в качестве управления использовать именно $u$.

Для системы (7) управление параметром $u$ возможно, т.к. для линеаризованной системы выполняется критерий Калмана.

Запишем функцию Гамильтона: \[ \tag{9} H = \lambda_1 (\mu(s) - D) x + \lambda_2 (- \mu(s) x - D s) + \lambda_2 D s_p + \lambda_2 D u. \] Так как $D \gt 0$, то оптимальное управление зависит от знака $\lambda_2$: \[ \tag{10} u = \begin{cases} - s_p, & \lambda_2 \lt 0, \\ \phantom{-} u_0, & \lambda_2 \gt 0. \end{cases} \]

Для исследования поведения $\lambda_2$ рассмотрим уравнения Эйлера-Лагранжа: \[ \tag{11} \begin{aligned} \dot \lambda_1 &= - \pd{H}{x} = - (\mu(s) - D) \lambda_1 + \lambda_2 \mu(s), \\ \dot \lambda_2 &= - \pd{H}{s} = - \pd{\mu(s)}{s} x \lambda_1 + \paren{ \pd{\mu(s)}{s} x + D } \lambda_2. \end{aligned} \] Если из второго вычтем первое, то получим \[ \dv{}{t} (\lambda_2 - \lambda_1) = \paren{D + \pd{\mu(s)}{s} x - \mu(s)} (\lambda_2 - \lambda_1). \] Пусть $\varphi = \lambda_2 - \lambda_1$, тогда \[ \dot \varphi(t) = \paren{ D + \pd{\mu(s)}{s} x - \mu(s) } \varphi. \] Решение этого уравнения: \[ \begin{aligned} \varphi(t) &= \varphi_0 e^{ \int\limits_{0}^{t} \paren{ D + \pd{\mu(s)}{s} x - \mu(s) } dt } = \\ &= \paren{ \lambda_2(0) - \lambda_1(0) } e^{ Dt + \int\limits_{0}^{t} \paren{ \pd{\mu(s)}{s} x - \mu(s) } dt }. \end{aligned} \] На промежутке от 0 до $t$ функция знак не меняет.

Выразим $\lambda_1$ как функцию от $\lambda_2$: \[ \lambda_1 = \lambda_1(\varphi, \lambda_2) = \lambda_2 - \varphi(t). \] Подставим это представление в уравнения (11). Второе уравнение имеет вид: \[ \dot \lambda_2 = D \lambda_2 + \pd{\mu(s)}{s} x \varphi(t). \] Решая это неоднородное уравнение по формуле Коши, получаем \[ \lambda_2(t) = e^{Dt} \paren{ \lambda_2(0) + \int\limits_{0}^{t} \pd{\mu(s)}{s} x \varphi(t) e^{-D \tau} d\tau }. \] Исследуем знак интеграла:

Итак, всё выражение под интегралом не меняет свой знак на промежутке $[0, t]$. Получается, что $\lambda_2(t)$ будет менять свой знак не более одного раза, следовательно, управление имеет не более одной точки переключения.

Посмотрим на то, как могут выглядеть оптимальные фазовые траектории. Пусть $L_1$ — семейство кривых (траекторий), соответствующих управлению $u = - s_p$ (то есть $s_0 = 0$). Тогда \[ \tag{12} \begin{aligned} \dot x &= \phantom{-} \paren{ \mu(s) - D } x, \\ \dot s &= - \mu(s) x - Ds. \end{aligned} \] Это случай, когда отсутствует подача субстрата в культиватор. Тогда \[ \tag{13} \dv{s}{x} = - \frac{\mu(s) x + D s}{\paren{ \mu(s) - D } x}. \] Так как \[ \mu(\tilde s) = D, \qquad \at{ \dv{s}{x} }_{s = \tilde s} = \infty, \] то $L_1$, проходя через кривую $s = \tilde s$, будут ей ортогональны.

Далее, \[ \mu(s) = \frac{s}{1 + s}, \] тогда \[ \frac{1}{\mu(s)} = 1 + \frac{1}{s}. \] Тогда $\mu(s)$ при увеличении $s$ монотонно возрастает. Следовательно, \[ s \gt \tilde s \implies \mu(s) \gt D \implies \dv{s}{x} \lt 0. \] И наоборот: \[ s \lt \tilde s \implies \mu(s) \lt D \implies \dv{s}{x} \gt 0. \]

Далее, складываем уравнение движения (12): \[ \dot x + \dot s = -D (x + s) \implies \dv{}{t} \paren{ x + s } = - D (x + s) \lt 0. \] Это означает, что \[ x + s \limto{t \to \infty} 0. \] Так как (из (12)) \[ \dv{s}{t} \lt 0, \implies s \limto{t \to \infty} 0, \] а поэтому \[ x \limto{t \to \infty} 0. \]

Далее, одна из параллельных прямых будет проходить через положение равновесия. Эта часть будет линией переключения. (что?) Так как есть только одна линия, проходящая через стационарную точку, то она и будет линией переключения. (проверь что тут написано)

Чтобы определить, в какую сторону движется точка по фазовым траекториям, заметим, что $\mu(s)$ монотонна, поэтому \[ \dot x \lt 0 \quad \text{при} \quad s \lt \tilde s. \] Получается, что проекция скорости точки на $x$ будет отрицательной. Значит, закон движения — по часовой стрелке.

Пусть теперь $u = u_0 = a - s_p$. Тогда \[ \tag{14} \begin{aligned} \dot x &= \phantom{-} \paren{ \mu(s) - D } x, \\ \dot s &= - \mu(s) x - Ds + D(s_p + u_0). \end{aligned} \] Проведя аналогичные исследования, выясняем, что траектории при таком управлении также ортогональны кривой $s = \tilde s$.

Если начинаем из $L_1$, то двигаемся по траектории до $P_2$, на ней переключаемся и уходим в стационарную точку. Аналогично с $L_2$ и $P_1$.

Выходит, что \[ \begin{aligned} \dot s &= - \mu(\tilde s) \tilde x - D \tilde s + D s_p + D u_0 = \\ &= - \cancel{D s_p} + \cancel{D \tilde s} - \cancel{D \tilde s} + \cancel{D s_p} + D u_0 = \\ &= D u_0 \gt 0. \end{aligned} \] Это означает, что $s$ возрастает при переходе через точку переключения.

Итак, $P_1 P_2$ — линия переключения, а управление имеет вид \[ u = \begin{cases} u_0, & (x,s) \in \Omega_+ \cup P_2, \\ - s_p, & (x,s) \in \Omega_- \cup P_1. \end{cases} \]

2023-11-09

Глава 4. Модели эпидемических процессов и управление ими

Математические модели эпидемий

  1. Механизмы заражения и их количественное описание

    Механизмы заражения:

    1. через контакт;
    2. через промежуточного хозяина;
    3. через среду.

    Считаем, что моделируем население. Оно делится на группы:

    1. восприимчивые;
    2. больные;
    3. иммунные.

    Предполагается, что численность населения не меняется, то есть надо учитывать даже умерших.

    Пусть $N$ — численность населения, $N_1, N_2, N_3$ — восприимчивые, заболевшие и иммунные соответственно.

    Рассмотрим модели для передачи.
    1. Контактная передача: \[ \Delta N_{i+1} = N_i N_j \Delta t. \] Пример: \[ \begin{aligned} \Delta N_1 &= - N_1 N_2 \Delta t, \\ \Delta N_2 &= N_1 N_2 \Delta t. \end{aligned} \]
    2. Переносчики (промежуточная популяция). Пусть $N_\text{п}$ — количество переносчиков, тогда \[ \begin{aligned} \Delta N_1 &= - N_1 N_\text{п} \Delta t, \\ \Delta N_2 &= N_1 N_\text{п} \Delta t, \\ \Delta N_\text{п} &= - N_1 N_\text{п} \Delta t. \end{aligned} \] В данном случае подразумевается, что переносчик умирает при контакте (например, всегда успеваем убить малярийного комара).
    3. Среда. Дополнительная характеристика — вирулентность $V$. Если она может изменяться, то \[ \Delta V = (c N_1 - m N_3) V \Delta t, \] где $c,m$ — эмпирически подобранные коэффициенты.
      Если вирулентность постоянна, то \[ \begin{aligned} \Delta N_1 &= - V N_1 \Delta t, & V &= \const, \\ \Delta N_2 &= \phantom{-} V N_2 \Delta t. \end{aligned} \]
  2. Модель без учёта изменчивости возбудителя (Кермака-МакКендрика)

    Представим её в разностном виде: \[ \begin{aligned} \Delta N_1 &= - V_{12} N_1 N_2 \Delta t, \\ \Delta N_2 &= V_{12} N_1 N_2 \Delta t - V_{23} N_2 \Delta t, \\ \Delta N_3 &= V_{23} N_2 \Delta t, \end{aligned} \] где $V_{12} = \const$ — вирулентность, $V_{23} = \const$ — некоторый весовой коэффициент.

    Переходя к пределам, получаем \[ \tag{1} \begin{aligned} \dot N_1 &= - k N_1 N_2, \\ \dot N_2 &= k N_1 N_2 - \beta N_2, \\ \dot N_3 &= \beta N_2. \end{aligned} \] Считаем, что $N = N_1 + N_2 + N_3 = \const$ и \[ \beta = \frac{1}{T}, \] где $T = \const$ — характерное время течения болезни — время, за которое количество заболевших уменьшается в $e$ раз.

    Начальные условия: \[ \begin{aligned} N_1(0) &= N_1^0 \gt 0, \\ N_2(0) &= N_2^0 \gt 0, \\ N_3(0) &= N_3^0 \geqslant 0. \end{aligned} \]

    Условие для начала эпидемии определяется вторым уравнением: $\dot N_2 \geqslant 0$, то есть \[ N_2 (k N_1 - \beta) \geqslant 0, \qquad N_1^0 \geqslant \frac{\beta}{k}. \]

    Построим приближённое решение. Введём обозначения: \[ x = N_1, \qquad y = N_2, \qquad z = N_3. \] Тогда $N = x + y + z$, а модель примет вид \[ \begin{aligned} \dot x &= - kxy, \\ \dot y &= kxy - \beta y, \\ \dot z &= \beta y, \end{aligned} \] начальные условия (примем $z_0 = 0$): \[ z(0) = z_0 =: 0, \quad x(0) = x_0 \gt 0. \] Разделим первое уравнение на третье, получим \[ \dv{x}{z} = - \frac{kx}{\beta} \implies \frac{dx}{x} = - \frac{k dz}{\beta}. \] Отсюда \[ \left. \ln x \right|_{x_0}^x = \left. - \frac{k}{\beta} z \right|_{z_0}^z. \] То есть \[ \ln \frac{x}{x_0} = - \frac{k}{\beta} z \implies x = x_0 e^{- \frac{kz}{\beta}}. \] Получили $x$ как функцию от $z$. Рассмотрим уравнение относительно $z$: \[ \dot z = \beta y. \] Так как \[ y = N - x - z, \] то \[ \begin{aligned} \dot z &= \beta (N - x - z) = \\ &= \beta \paren{ N - z - x_0 e^{- \frac{kz}{\beta}} }. \end{aligned} \] Чтобы получить приближённое значение, разложим экспоненту по степеням (до второй степени): \[ \begin{aligned} \dot z &= \beta N - \beta x_0 \paren{ 1 - \frac{k}{\beta} z + \frac{k^2}{2 \beta^2} z^2 } - z \beta = \\ &= \beta(N - x_0) + (k x_0 - \beta) z - x_0 \frac{k^2}{2 \beta} z^2. \end{aligned} \] Получили уравнение Риккати с постоянными коэффициентами. Решая его, получаем \[ z(t) = \frac{\beta^2}{k^2 x_0} \left[ \frac{k}{\beta} x_0 - 1 + \sqrt{q} \th \paren{ \frac{\sqrt{q}}{2} \beta t - \varphi } \right], \] где \[ \begin{aligned} \th x &= \frac{e^x - e^{-x}}{e^x + e^{-x}}, \\ \varphi &= \th^{-1} \paren{ \frac{k x_0 / \beta - 1}{\sqrt{q}} }, \\ \sqrt{q} &= \sqrt{ \paren{ \frac{k x_0}{\beta} - 1 }^2 + 2 x_0 y_0 \frac{k^2}{\beta^2} }, \end{aligned} \] причём $y_0 = N - x_0 - z_0$.

  3. Модель с учётом изменчивости возбудителя
    Рассматриваем модель \[ \tag{2} \begin{aligned} \dot N_1 &= - k V N_1 N_2, \\ \dot N_2 &= k V N_1 N_2 - \beta N_2, \\ \dot N_3 &= \beta N_2, \\ \dot V &= (c N_1 - m N_3) V. \end{aligned} \] Также $N = N_1 + N_2 + N_3 = \const$ и \[ \begin{aligned} N_1(0) &\gt 0, \\ N_2(0) &\gt 0, \\ N_3(0) &\geqslant 0, \\ V(0) &\gt 0. \end{aligned} \]

    В случае дрейфа можно учесть падение коллективного иммунитета: \[ \begin{aligned} \dot N_1 &= - k V N_1 N_2 + \frac{N_3}{T_D}, \\ \dot N_2 &= k V N_1 N_2 - \beta N_2, \\ \dot N_3 &= \beta N_2 - \frac{N_3}{T_D}, \\ \dot V &= (c N_1 - m N_3) V. \end{aligned} \] где $T_D$ — характерное время утраты коллективного иммунитета.

    Если взять за $T_I$ характерное время длительности иммунитета, то можно положить \[ q = \frac{1}{T_I} + \frac{1}{T_D} = \frac{1}{T_r} \implies T_r = \frac{1}{q}. \]

  4. Учёт латентной фазы (модель Ренальда)

    Введём дополнительную группу $N_L$ — количество латентно заражённых, то есть \[ N = N_1 + N_2 + N_3 + N_L. \] Тогда получаем следующую систему: \[ \tag{3} \begin{aligned} \dot N_1 &= - k N_1 N_2 + \frac{N_3}{T_r}, \\ \dot N_L &= \phantom{-} k N_1 N_2 - \frac{N_L}{T_L}, \\ \dot N_2 &= - N_2 \beta + \frac{N_L}{T_L}, \\ \dot N_3 &= \phantom{-} N_2 \beta - \frac{N_3}{T_r}, \end{aligned} \] где $T_L$ — характерная длительность инкубационного периода. Предполагаем, что латентно заражённые не являются заразными.

    В данном случае условие начала эпидемии будет \[ \dot N_2 \geqslant 0 \implies N_L^0 \geqslant T_L \beta N_2^0. \]

Модели управления в паразитарных системах

  1. Способы управление эпидемическим процессом

    Инфекции делятся на управляемые и неуправляемые.

    Возможные типы управления:
    1. вакцинация;
    2. изоляция;
    3. профилактика.
  2. Модели с изоляцией
    1. Карантин.
      Управление $u_k$ — интенсивность изоляции (количество изолированных в сутки), $N_k$ — число изолированных. Считаем $V = \const$. Тогда \[ \begin{aligned} \dot N_1 &= -k V N_1 N_2 - u_k, \\ \dot N_2 &= k V N_1 N_2 - \beta N_2, \\ \dot N_3 &= \beta N_2, \\ \dot N_k &= u_k, \end{aligned} \] где $N = N_1 + N_2 + N_3 + N_k$.
    2. Изоляция больных.
      Управление $u_b$ — интенсивность изоляции. Оно не будет превышать количество заболевших в сутки. Тогда $N_2'$ — количество изолированных больных. Получаем систему \[ \begin{aligned} \dot N_1 &= -k V N_1 N_2, \\ \dot N_2 &= k V N_1 N_2 - \beta N_2 - u_b, \\ \dot N_3 &= \beta N_2 + \beta N_2', \\ \dot N_2' &= u_b - \beta N_2'. \end{aligned} \]
  3. Модели с профилактикой и дезинфекцией
    1. Химия.
      Введём доли тех, кто пользуется и кто не пользуется профилактикой: \[ K + K' = 1, \] где $K$ — не пользуются, а $K'$ — пользуются. Обозначим \[ \gamma = \frac{1}{T_L}, \qquad \beta = \frac{1}{T}, \] тогда \[ \begin{aligned} \dot N_1 &= - k V N_1 N_2 + \gamma K' N_L, \\ \dot N_L &= k V N_1 N_2 - \gamma (K + K') N_L, \\ \dot N_2 &= \gamma K N_L - \beta N_2, \\ \dot N_3 &= \beta N_2. \end{aligned} \]
    2. Дезинфекция.
      Будем уменьшать коэффициент $k$.
    3. Интерферон.
      Пусть \[ \dot V = (c N_1 - m N_2) V. \] Тогда интерферонопрофилактика обрушивает положительную обратную связь, то есть $c \to \min$. Если же $V = \const$, то $k \to \min$. Обычно $\min = 0$.
  4. Модели с вакцинацией
    1. Простая вакцинация.
      Управление — количество вакцинированных в сутки. Вакцинация может использоваться только до момента начала эпидемии. \[ u_v(t) = \begin{cases} \tilde u_v \gt 0, & t \in [0, T_v], \\ 0, & t \gt T_v, \end{cases} \] где $T_v$ — время начала эпидемии.

      Модель выглядит следующим образом: \[ \tag{4} \begin{aligned} \dot N_1 &= - k V N_1 N_2 - u_v, \\ \dot N_2 &= k V N_1 N_2 - \beta N_2, \\ \dot N_3 &= \beta N_2 + u_v, \end{aligned} \] причём $V = \const$ и $N = N_1 + N_2 + N_3 = \const$.

    2. Комбинированная вакцинация.
      Рассмотрим, например, грипп. Пусть вакцинируемся сразу и от гриппа A, и от гриппа B. Эпидемию считаем кратковременной, то есть человек не успевает заразиться другим видом.

      Пусть $N_{3A}, N_{3B}$ — число переболевших гриппом $A$ и $B$, а $\delta_A, \delta_B$ — доли тех, кто защитился от $A$ и $B$, причём \[ \delta_A + \delta_B = 1. \] Тогда \[ \begin{aligned} \dot N_1 &= - k V N_1 N_2 - u_{AB}, \\ \dot N_2 &= k V N_1 N_2 - \beta_A N_2 - \beta_B N_2, \\ \dot N_{3A} &= \beta_A N_2 + \delta_A u_{AB}, \\ \dot N_{3B} &= \beta_B N_2 + \delta_B u_{AB}, \end{aligned} \] причём $N = N_1 + N_2 + N_{3A} + N_{3B} = \const$.

2023-11-16

Оптимизация вакцинации на основе принципа максимума

Рассмотрим модель \[ \tag{5} \begin{aligned} \dot N_1 &= - K V N_1 N_2 - u, \\ \dot N_2 &= K V N_1 N_2 - \beta N_2, \\ \dot N_3 &= \beta N_2 + u, \\ \dot V &= (c N_1 - m N_3) V, \end{aligned} \] причём $N = N_1 + N_2 + N_3 = \const$. Начальные условия: \[ \begin{aligned} N_1(0) &= N_1^0 \gt 0, \\ N_2(0) &= N_2^0 \gt 0, \\ N_3(0) &= N_3^0 \geqslant 0, \\ V(0) &= V_0 \gt 0. \end{aligned} \] Проводим замену: \[ N_1 = N - N_2 - N_3, \] тогда \[ \tag{6} \begin{aligned} \dot N_2 &= K V N_2 \paren{ N - N_2 - N_3 } - \beta N_2, \\ \dot N_3 &= \beta N_2 + u, \\ \dot V &= \paren{ a - c N_2 - b N_3 } V, \end{aligned} \] где \[ a = c N, \qquad b = c + m. \] Критерий качества — экономия средств, потраченных на вакцинацию и лечение. Пусть $C_L$ — цена лечения, $C_V$ — цена вакцинации, причём \[ C_L \gg C_V. \] Пусть $t = 0$ — момент начала вакцинации, $T'$ — момент окончания вакцинации, $T$ — момент окончания эпидемии. Тогда критерий качества принимает вид \[ \tag{7} J = \int\limits_{0}^{T} (C_V u + C_L V K N_2 (N - N_2 - N_3)) dt \to \min. \]

Ограничение на управление: \[ \tag{8} U = \begin{cases} u_{\min} \leqslant u \leqslant u_{\max}, & t \in [0, T'), \\ u = 0, & t \in [T', T], \end{cases} \] где \[ u_{\max} = \frac{\delta N}{T_0}, \qquad \delta = 0.75, \] а $T_0$ — средняя длительность эпидемии по многолетним наблюдениям.

Составляем функцию Гамильтона: \[ \begin{aligned} H &= - C_V u - C_L V K N_2 (N - N_2 - N_3) + \\ &\phantom{=} + \lambda_1 \left[ K V N_2 (N - N_2 - N_3) - \beta N_2 \right] + \\ &\phantom{=} + \lambda_2 \left[ \beta N_2 + u \right] + \lambda_3 \paren{ a - c N_2 - b N_3 } V. \end{aligned} \] Отсюда следует, что \[ \max H = \max (\lambda_2 u - C_V u) = \max (\lambda_2 - C_V) u. \] Понятно, что \[ \tag{9} u_o = \begin{cases} u_{\min}, & \lambda_2 - C_V \lt 0, \\ u_{\max}, & \lambda_2 - C_V \gt 0. \end{cases} \] Для нахождения $\lambda_2$ выпишем уравнения Эйлера-Лагранжа: \[ \tag{10} \begin{aligned} \dot \lambda_1 &= - \pd{H}{N_2}, \\ \dot \lambda_2 &= - \pd{H}{N_3}, \\ \dot \lambda_3 &= - \pd{H}{V}. \end{aligned} \] Воспользуемся условиями трансверсальности: \[ \left. \paren{\lambda \Delta x - H \Delta t} \right|_{0}^{T} = 0. \] Считаем, что $t = 0$ и $T$ — фиксированы, поэтому $\left. H \Delta t \right|_{0}^{T} = 0$. Заметим, что $x_0$ фиксирован, поэтому $\at{\Delta x_0}_{t = 0} = 0$. Отсюда \[ \lambda(T) \Delta x(T) = 0. \] Так как $\Delta x(T)$ — независимые вариации, то для выполнения этого равенства требуем \[ \lambda(T) = 0. \]

Решать будем численным алгоритмом Черноусько-Крылова. Пусть $u_0(t) = u^0$ — некоторое начальное приближение, которое выбирается на основе многолетнего опыта.

  1. Решаем в прямом интегрировании (6) и в обратном интегрировании (10) при $u = u^0$. Получаем значения для оптимальной траектории (из (6)) \[ x^* = x^*(N_2^*, N_3^*, V^*), \] а с помощью пары $(x^*, u^0)$ получаем множители Лагранжа (из (10)): \[ \lambda^* = \lambda^*(t). \]
  2. Из функции Гамильтона определяем точки переключения: \[ \pd{H}{u} = 0, \implies \lambda_2 = C_V. \] Понятно, что $\lambda^* = (\lambda_1^*, \lambda_2^*, \lambda_3^*)$. Используя $\lambda_2^*(t) = C_V$, находим новые точки переключения для $u^1(t)$. Проверяем, насколько новое приближение $u^1(t)$ улучшает минимизацию функционала (7).
  3. Возвращаемся к шагу 1, пока либо значения перестали сильно изменяться, либо значения критерия качества изменяются в пределах заранее заданной погрешности.

Глава 5. Математическая теория управляемых экономических систем

Пункт 1. Примеры статических и динамических моделей управляемых экономических систем

Будем рассматривать функцию полезности и производственную функцию.

Функция полезности описывает потребление. По сути представляет собой набор предпочтений того или иного потребителя.

  1. Оптимизация экономической статики
    Рассмотрим набор $x = (x_1, x_2, \dots, x_n)$.
    Будем говорить, что набор $x$ предпочтительнее набора $x'$, если \[ f(x) \gt f(x'), \] где $f(x)$ — функция полезности.
    Свойства функции полезности:
    1. Если $\tilde x_i \gt x_i$, то \[ f(x_1, \dots, \tilde x_i, \dots, x_n) \gt f(x_1, \dots, x_i, \dots, x_n). \] Возрастание потребления одного продукта из набора влечёт за собой возрастание функции полезности. То есть \[ \pd{f}{x_i} \gt 0, \qquad i = \overline{1, n}. \]
      Частные производные \[ \pd{f}{x_i}, \qquad i = \overline{1, n}. \] называются предельными полезностями.
    2. Предельная полезность каждого продукта уменьшается с увеличением объёма потребления этого продукта: \[ \pdv2{f}{x_i} \lt 0. \]
    3. Предельная полезность каждого продукта увеличивается, если растёт объём потребления другого продукта: \[ \ppdv{f}{x_i}{x_j} \gt 0. \]
    4. Линии равного уровня полезности \[ f(x) = \const \] называются линиями безразличия.
    Задача максимизации полезности \[ f(x) \to \max \] при бюджетном ограничении \[ \sum\limits_{i=1}^{n} x_i \leqslant I, \qquad x_i \gt 0. \]
    Набор $x_0$, являющийся решением поставленной оптимизационной задачи, называют локальным рыночным равновесием.
    Рассмотрим набор из двух благ: $x_1, x_2$. В этом случае функция полезности может быть записана в виде \[ f(x_1, x_2) = a_1 \lg(x_1 - \overline x_1) - a_2 \lg(x_2 - \overline x_2), \] где $\overline x_1, \overline x_2$ — минимальные средства, необходимые для приобретения блага.

    Рассмотрим оптимизацию производственной функции.

    Производственная функция: \[ y = f(K, L), \] где $y$ — количество произведённого товара, $K$ — капитал, а $L$ — затраты на труд.
    5 рабочих за 1 час на 3 станках изготавливают 20 деталей, тогда
    • $y = 20$ товаров,
    • $L = 5$ человекочасов,
    • $K = 3$ станкочаса.
    Производственная функция Кобба-Дугласа: \[ y(K, L) = a_0 K^{a_1} L^{a_2}, \qquad K, L \gt 0, \] где $a_0, a_1, a_2 = \const$ — определяются эмпирическим образом.
    Производственная функция Леонтьева: \[ y(K, L) = \min \paren{ \frac{K}{a}, \frac{L}{b} }, \] где $a, b = \const$ — определяются эмпирическим образом.
    Свойства производственной функции:
    1. $f(x_1, \dots, 0, \dots, x_n) = 0$.
    2. Увеличение любых затрат влечёт увеличение производства: \[ \pd{f}{x_i} \gt 0. \]
    3. С увеличением затрат на ресурс уменьшается темп увеличения производства: \[ \pdv2{f}{x_i} \lt 0. \]
    4. Эффективность не падает при увеличении затрат на что-нибудь другое: \[ \ppdv{f}{x_i}{x_j} \geqslant 0. \]
    Рассмотрим коэффициент научного прогресса. Например, для функции Кобба-Дугласа \[ y = a_0 K^{a_1} L^{a_2} e^{\lambda t}, \qquad \lambda \gt 0. \] Например, для позднего СССР $\lambda \approx 0.03$.
    Примеры задач в статике.
    Пусть $c_k$ — цена единицы капитализации, а $c_L$ — цена труда, то бюджетное ограничение представимо в виде \[ c_k K + c_L L \leqslant G. \] Тогда надо найти максимум \[ y(K, L) \to \max \] при соблюдении ограничений \[ c_k K + c_L L \leqslant G. \]
    Пусть $c_0$ — цена за единицу выпускаемого товара, тогда можно поставить задачу максимизации прибыли: \[ F(K, L) = c_0 y(K, L) - c_K K - c_L L \to \max \] при условии \[ c_k K + c_L L \leqslant G. \]
    Рассмотрим $n$ предприятий, у каждого своя производственная функция: \[ y_j = f_j(x_1^j, \dots, x_m^j), \qquad j = \overline{1,n}, \] причём ресурсы ограничены: \[ \sum\limits_{j=1}^{m} x_i^j \leqslant X_i, \qquad i = \overline{1,n}. \] Тогда можно поставить задачу максимизации общего производства: \[ \sum\limits_{j=1}^{n} y_j \to \max. \]
  2. 2023-11-23

  3. Оптимизация экономической динамики
    1. Динамика несвязанных секторов экономики
      Предполагается, что экономика состоит из нескольких несвязанных секторов. Каждый из секторов имеет свою производственную функцию. Как они задаются — не конкретизируем.

      Введём коэффициенты: $\alpha_i$ — обозначают ту долю произведённого продукта, идущую на накопление капитала, и $\beta_i$ — выбывающую из накопления долю капитала.

      Тогда \[ \begin{aligned} \dot K_1 &= \alpha_1 y_1 - \beta_1 K_1, \\ \dot K_2 &= \alpha_2 y_2 - \beta_2 K_2, \\ \dot K_3 &= \alpha_3 y_3 - \beta_3 K_3. \end{aligned} \] Начальные условия: \[ K_i(0) = K_i^0. \] Получаем автономные уравнения, но их совместное решение даёт картину всей экономики.

      Находим равновесное положение и решаем уравнения по отдельности. Часто предполагают, что $L_i = \const$.

    2. Динамика связанных секторов экономики
      Рассмотрим \[ \tag{1} \dot K(t) = \alpha y(K, L) - \beta K, \qquad K(t) \in \R^3. \] Причём \[ \alpha = \begin{pmatrix} \alpha_{11} & \alpha_{12} & \alpha_{13} \\ \alpha_{21} & \alpha_{22} & \alpha_{23} \\ \alpha_{31} & \alpha_{32} & \alpha_{33} \end{pmatrix}, \] где $\alpha_{ij}$ — перетекание доходов сектора $j$ в сектор $i$.

Пункт 2. Примеры постановок задач оптимального управления для системы $(1)$

  1. Задача с фиксированными концами (увеличение капитала до заданного значения)

    Пусть \[ K(0) = K^0, \qquad K(T) = K^T. \] Время $T$ считается фиксированным. Рассматриваем уравнение (1) \[ \dot K(t) = \alpha y(K, L) - \beta K. \] Пусть $K, L \in \R$. В качестве управления рассматриваем затраты на труд: $u = L$, причём они ограничены: \[ 0 \leqslant L \leqslant L_{\max}. \] В качестве критерия оптимальности рассматриваем также затраты на труд: \[ J = \int\limits_{0}^{T} L^2(t) dt \to \min. \]

    Решаем задачу принципом максимума. Составляем функцию Гамильтона \[ H = H(K, L, \lambda) = -L^2 + \lambda \paren{ \alpha y(K, L) - \beta K }. \] Рассмотрим уравнение Эйлера-Лагранжа \[ \dot \lambda = -\pd{H}{K}. \] Понятно, что они зависят от вида производственной функции $y(K, L)$.

  2. Задача со свободным правым концом (максимизация потребления)

    Фиксируем левый конец $K(0) = K^0$ и время $T = \const$. Конечное значение для капитала не фиксировано. В качестве управления используем коэффициент капитализации $u = \alpha$. Тогда критерий качества имеет следующий вид: \[ J = \int\limits_{0}^{T} (1 - \alpha) y(K, L) dt \to \max. \]

    Ограничения на управление: \[ \alpha_{\min} \leqslant \alpha \leqslant \alpha_{\max}. \] Также будем предполагать, что затраты на труд постоянны: $L = \const$.

    Гамильтониан имеет следующий вид: \[ H = y(K, L) - \lambda \beta K + \alpha \paren{ \lambda y(K, L) - y(K, L) }. \] Он линеен по управлению, поэтому максимальное значение гамильтониана достигается на концах огранчинеий, то есть при $\alpha_{\min}$ или при $\alpha_{\max}$.

    Общее условие трансверсальности: \[ \at{ \paren{ \lambda \Delta K - H \Delta t } }_{0}^T = 0. \] Оно принимает вид \[ \lambda(T) \Delta K(T) = 0. \] Из независимости вариации следует, что \[ \lambda(T) = 0. \]

  3. Задача с подвижным правым концом

    Рассмотрим задачу оптимального управления для двух секторов экономики: \[ \begin{aligned} \dot K_1 &= \alpha_{11} y_1(K_1, L_1) + \alpha_{12} y_2(K_2, L_2) - \beta_1 K_1, \\ \dot K_2 &= \alpha_{21} y_1(K_1, L_1) + \alpha_{22} y_2(K_2, L_2) - \beta_2 K_2. \end{aligned} \] Начальный капитал: \[ K_1(0) = K_1^0, \qquad K_2(0) = K_2^0. \] Производственные функции: \[ \begin{aligned} y_1 &= a_0^1 K_1^{a_1^1} L_1^{a_2^1}, \\ y_2 &= a_0^2 K_2^{a_1^2} L_1^{a_2^2}. \end{aligned} \] В качестве управления рассматриваем вектор трудозатрат \[ u = L = (L_1, L_2). \] Критерий оптимальности: \[ J = \int\limits_{0}^{T} (L_1^2(t) + L_2^2(t)) dt \to \min. \] Время $T = \const$ — фиксировано.

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

    Гамильтониан: \[ \begin{aligned} H &= -L_1^2 - L_2^2 + \lambda_1 \paren{ \alpha_{11} y_1(K_1, L_1) + \alpha_{12} y_2(K_2, L_2) - \beta_1 K_1 } + \\ &\phantom{= -L_1^2 - L_2^2} + \lambda_2 \paren{ \alpha_{21} y_1(K_1, L_1) + \alpha_{22} y_2(K_2, L_2) - \beta_2 K_2 }. \end{aligned} \]

    Необходимые условия оптимальности: \[ \pd{H}{L_1} = 0, \qquad \pd{H}{L_2} = 0. \]

    Можно найти оптимальное управление: \[ \begin{aligned} L_1^0 &= \left[ \frac{1}{2} a_0^1 a_2^1 K_1^{a_1^1} \paren{ \lambda_1 \alpha_{11} + \lambda_2 \alpha_{21} } \right]^{a_2^1 - 2}, \\ L_2^0 &= \left[ \frac{1}{2} a_0^1 a_2^2 K_2^{a_1^2} \paren{ \lambda_1 \alpha_{12} + \lambda_2 \alpha_{22} } \right]^{a_2^2 - 2}. \end{aligned} \]

    Получили оптимальное управление как функцию от $K$ и $\lambda_1, \lambda_2$. Для их решения надо решать систему уравнений: \[ \dot \lambda_i = - \pd{H}{K_i}. \] Проблемы возникают, тк нужны граничные значения для $\lambda_1, \lambda_2$.

    Если, например, в конечный момент времени должно выполняться соотношение \[ K_1 = \varphi(K_2), \] то выписываем условия трансверсальности: \[ \lambda(T) \Delta K(T) = 0. \] Так как \[ \dv{K_1}{K_2} = \varphi'(K_2), \] то, переходя к вариациям, получаем зависимость \[ \Delta K_1 = \Delta K_2 \varphi'(K_2). \] Подставляем это в условие трансверсальности: \[ \lambda_1(T) \Delta K_1(T) + \lambda_2(T) \Delta K_1(T) = 0. \] Получаем уравнение относительно одной вариации, откуда \[ \lambda_1(T) \varphi'(K_2 (T)) + \lambda_2(T) = 0. \] Находим соотношение между $\lambda_1(T), \lambda_2(T)$, благодаря которому можно решить уравнения Эйлера-Лагранжа.

Глава 6. Оптимальное импульсное управление

§1. Вводные понятия

1. Расширение классической задачи оптимизации
Если в момент времени $\theta$ происходит единичный импульс, то обозначать это будем через дельта-функцию Дирака $\delta( t - \theta )$.

Будем использовать следующий вид управления: \[ \tag{1} \nu(t) = u(t) + \sum\limits_{i=0}^{k} c_i \delta(t - \theta_i), \] где $u(t) \in PC$, а второе слагаемое заключает в себе возможные импульсы в моменты времени $\theta_i$ величиной $c_i$. Предполагаем, что \[ 0 \leqslant \theta_0 \lt \dots \lt \theta_k \leqslant T. \] Тогда для поиска оптимального управления можем оптимизировать:

Считаем, что $\nu \in U$, причём $U$ неограниченное. Будем рассматривать случаи $U = \R$ и $U = \R_+$.

Как это управление воздействует на правые части системы? Оно может вызывать разрывы.

2. Структура разрывных траекторий
  1. Введём множество $D$ — множество вектор-функций $x(t) \in D$, являющихся решениями задачи оптимального управления, которые на промежутке времени $t \in [0; T]$ удовлетворяют условиям:
    1. имеют на этом промежутке не более конечного числа разрывов первого рода, причём разрывы не исключены в точках $t = 0$ и $t = T$;
    2. в промежутках между разрывами $x(t)$ непрерывна и кусочно-дифференцируема, то есть имеет свойства обычной траектории;
    3. для всех точек разрыва $x(t)$ непрерывна справа при $(0, T]$, то есть \[ x(\theta) = x(\theta_+), \qquad x(T) = x(T_+), \] а при $t = 0$ \[ x(0) = x(0_-). \]
  2. Пусть задана управляемая система \[ \tag{2} \dot x(t) = f(t, x(t)) + g(t, x(t)) \nu, \qquad x(0) = x^0, \] причём $f, g$ — гладкие и $x(t) \in D$. Эта система линейна по управлению.

    Введём вспомогательную систему: \[ \tag{3} \dv{z_i}{\tau} = g_i(\theta_j, z_i) c_j, \qquad z_i(0) = x_i(\theta_{j-}), \qquad \tau \in [0, 1]. \]

    Систему (3) называют предельной.

    Она вводится для того, чтобы убедиться, что это действительно импульс, а не просто разрыв.

    Условие допустимости скачка: \[ \tag{4} z_i(1) = x_i(\theta_{j+}). \] Если оно выполнено, то действительно в момент времени $t = \theta_j$ происходит импульс, где $j = \overline{1,k}$.
  3. Функцию $x \in D$ называют обобщённой траекторией системы (2), соответствующей импульсному управлению (1), если в любой момент импульса $\theta_i$ функция $x(t)$ удовлетворяет условию допустимости скачка (4), а на интервалах между импульсами — обычной дифференциальной системе \[ \dot x(t) = f(t, x(t)) + g(t, x(t)) u(t). \]
Задача максимизации дисконтированного суммарного потребления модели роста капитала Рамсея.
Дисконтирование — процесс оценивания будущих доходов в настоящем.
Рассмотрим уравнение баланса: \[ Y(t) = F(K(t)) = C(t) + I(t), \] где Уравнение развития капитала: \[ \dot K(t) = I(t), \qquad K(0) = K^0. \] Критерий качества: \[ J = \int\limits_{0}^{T} e^{-r t} C(t) dt \to \max, \] где $r$ — ставка дисконта. Момент времени $T = \const$ — фиксирован.

В качестве управления — $I(t)$. Инвестиции ограничены сверху, но также допускается неограниченное изъятие, то есть $I(t)$ может быть отрицательным.

На $[0, T]$ должно быть выполнено условие (фазовое ограничение): \[ K(t) \geqslant 0. \]

Преобразуем эту модель. Преобразуем критерий Лагранжа в критерий Майера путём ввода дополнительной переменной, а также перейдём к другому управлению.

С помощью уравнения баланса избавимся от управления: \[ I(t) = F(K(t)) - C(t). \] Введём дополнительную переменную: \[ P(t) = \int\limits_{0}^{t} e^{-r s} C(s) ds. \] Тогда критерий качества превратится в \[ J = P(T) \to \max. \] Далее, с учётом уравнения баланса, можно переписать систему в следующем виде: \[ \dot P(t) = e^{-rt} C(t), \qquad P(0) = 0. \] Уравнение развития капитала: \[ \dot K(t) = F(K(t)) - C(t), \qquad K(0) = K^0. \] Фазовые ограничения: \[ K(t) \geqslant 0, \qquad C(t) \geqslant 0. \] В качестве управления рассматриваем $C(t)$.

Пусть в момент времени $\theta$ имеется разрыв величины $n$: \[ n \delta(t - \theta). \] Составим две предельные системы: для $P(t)$ и для $K(t)$: \[ \begin{aligned} \dv{z_P}{\tau} &= e^{-r \theta} n, & z_P(0) &= P(\theta_-), \\ \dv{z_K}{\tau} &= -n, & z_K(0) &= K(\theta_-). \end{aligned} \]

Интегрируем эти уравнения: \[ \begin{aligned} z_P(\tau) &= P(\theta_-) + e^{-r \theta} n \tau, \\ z_K(\tau) &= K(\theta_-) - n \tau. \end{aligned} \] Правые пределы траекторий должны равняться значениям траекторий в момент времени $\tau = 1$, то есть \[ \begin{aligned} P(\theta_+) &= z_P(1) = P(\theta_-) + e^{-r\theta} n, \\ K(\theta_+) &= z_K(1) = K(\theta_-) - n. \end{aligned} \] Отсюда следует, что \[ \begin{aligned} P(\theta_+) - P(\theta_-) &= e^{-r\theta} n, \\ K(\theta_+) - K(\theta_-) &= n. \end{aligned} \]

2023-11-30

§2. Принцип максимума в задачах импульсного управления

1. Постановка задачи оптимизации

Будем рассматривать управляемую систему \[ \tag{1} \dot x = f(x(t), t) + g(x(t), t) \nu, \] где \[ \tag{2} \nu(t) = u(t) + \sum\limits_{i=1}^{k} c_i \delta(t - \theta_i), \qquad 0 \leqslant \theta_0 \lt \dots \lt \theta_k \leqslant T. \] Будем считать, что задано начальное условие \[ x(0) = x^0. \] Ограничение на управление: \[ \tag{3} \nu \in U, \qquad c_i \in U, \] где $U = \mathbb{R}$ или $U = \mathbb{R}_+$.

Условие допустимости скачка в момент времени $t = \theta_j$: \[ \tag{4} \dv{z_i}{\tau} = g_i(\theta_j, z_i) c_j, \qquad \tau \in [0; 1], \] где \[ z_i(0) = x_i(\theta_{j-}), \qquad z_i(1) = x_i(\theta_{j+}). \]

Рассмотрим целевой функционал в форме Больца: \[ \tag{5} \int\limits_{0}^{T} F_0(t, x(t), \nu(t)) dt + \varphi_0(x(T)) \to \min, \] где \[ \tag{6} F_0 = f(t, x(t)) + g_0(t, x(t)) \nu(t). \]

2. Гамильтониан и уравнения Эйлера-Лагранжа

Рассмотрим гамильтониан \[ \tag{7} H = -F_0 + \lambda \dot x. \] Тогда уравнения Эйлера-Лагранжа записываются так: \[ \tag{8} \dot \lambda = - \pd{H(t, x(t), \lambda(t), \nu(t))}{x}, \] а условия трансверсальности — так: \[ \tag{9} \lambda(T) = - \pd{\varphi_0(x(T))}{x}. \] Заметим, что \[ \lambda(T) = \lambda(T^+), \qquad x(T) = x(T^+). \]

Составим предельную систему (4) для проверки допустимости скачка в момент времени $t = \theta_j$. В этой системе справедливы следующие соответствия переменных: \[ x \to z, \quad t \to \theta, \quad \nu \to c, \quad \lambda \to p, \quad H \to h. \] Пусть \[ h(t, x(t), \lambda(t), \nu(t)) = \lambda g(t, x) \nu - g_0(t, x) \nu, \] тогда \[ \dv{p_i}{\tau} = - \pd{h_i(\theta_j, z_i, p_i, c_j)}{z_i}, \qquad \tau \in [0,1], \] где $z_i$ — решение системы (4). Начальное условие: \[ p_i(0) = \lambda_i(\theta_{j-}). \] Условие допустимости скачка: \[ p_i(1) = \lambda_i(\theta_{j+}). \]

Пусть задана система уравнений движения: \[ \begin{aligned} \dot x_1 &= x_2 + \nu, \\ \dot x_2 &= - \nu, \end{aligned} \] начальные условия: \[ x_1(0) = x_{10}, \qquad x_2(0) = x_{20}, \] ограничение на управление: $U = \mathbb{R}$, целевой функционал: \[ J = 0.5 \int\limits_{0}^{T} x_1^2 dt. \]

Составим гамильтониан: \[ H = - 0.5 x_1^2 + \lambda_1 (x_2 + \nu) - \lambda_2 \nu \] и уравнения Эйлера-Лагранжа: \[ \begin{aligned} \dot \lambda_1 &= x_1, \\ \dot \lambda_2 &= - \lambda_1. \end{aligned} \] Заметим, что уравнения Эйлера-Лагранжа не зависят от $\nu$, поэтому предельная система для $\lambda$ не нужна — имульса не будет.

Составим предельную систему в момент времени $t = \theta$: \[ \begin{aligned} \dv{z_1}{\tau} &= \phantom{-} c, \\ \dv{z_2}{\tau} &= -c. \end{aligned} \] После интегрирования получаем \[ \begin{aligned} z_1(\tau) &= \phantom{-} c \tau + z_1(0), \\ z_2(\tau) &= -c \tau + z_2(0), \end{aligned} \] но из начальных условий известно, что \[ z_1(0) = x_1(\theta_{-}), \qquad z_2(0) = x_2(\theta_{-}), \] поэтому \[ \begin{aligned} z_1(\tau) &= \phantom{-} c \tau + x_1(\theta_{-}), \\ z_2(\tau) &= -c \tau + x_2(\theta_{-}). \end{aligned} \] Так как $\tau \in [0, 1]$, то рассмотрим значение решения предельной системы в момент времени $\tau = 1$ \[ \begin{aligned} z_1(1) &= \phantom{-} c + x_1(\theta_{-}), \\ z_2(1) &= -c + x_2(\theta_{-}). \end{aligned} \] Из условий допустимости скачка \[ z_1(1) = x_1(\theta_{+}), \qquad z_2(1) = x_2(\theta_{+}) \] следует, что \[ \begin{aligned} x_1(\theta_{+}) &= \phantom{-} c + x_1(\theta_{-}), \\ x_2(\theta_{+}) &= -c + x_2(\theta_{-}), \end{aligned} \implies \begin{aligned} x_1(\theta_{+}) - x_1(\theta_{-}) &= \phantom{-} c, \\ x_2(\theta_{+}) - x_2(\theta_{-}) &= -c. \end{aligned} \] Значит, в момент времени $t = \theta$ допустим скачок по $x_1$ на $c$, а по $x_2$ — на $-c$.

2023-12-07

3. Принцип максимума
Пусть оптимальное управление имеет вид \[ \nu^* = u^*(t) + \sum\limits_{i=1}^{k} c_i^* \delta (t - \theta_i^*). \] Тогда $x^*(t), \lambda^*(t)$ — решение основной и сопряжённой задачи.
Функция переключения: \[ H_\nu = \pd{H}{\nu}. \]
Необходимое условие оптимальности (ПМП) для задачи (1)-(5): должны выполняться условия:
  1. Для любого $t \in [0, T]$, за исключением точек разрыва обычной части траектории и за исключением моментов импульсов, обычная часть управления $u^*(t)$ максимизирует следующую линейную функцию: \[ H_\nu(t, x^*(t), \lambda^*(t)) \cdot u^*(t) = \max_{u \in U} H_\nu(t, x^*(t), \lambda^*(t)) u(t). \]
  2. В каждой точке импульса $\theta_i^*$ функция переключения равна нулю (в силу линейности): \[ H_\nu = H_\nu(\theta_i^*, x^*(\theta_i^*), \lambda^*(\theta_i^*)) = 0. \]
  3. Функция переключения $H_\nu(t, x^*(t), \lambda^*(t))$ непрерывна на промежутке $[0, T]$ и дифференцируема всюду, кроме точек приложения импульсов $\theta_i^*$.
  4. Функция \[ H_0(t, x(t), \lambda(t)) = \lambda f(t, x) - f_0(t, x) \] непрерывна на открытом промежутке $(0, T)$ и дифференцируема на нём всюду, кроме точек импульсов, причём \[ \pd{H_0(t, x^*, \lambda^*)}{t} = \dv{H_0(t, x^*, \lambda^*)}{t}. \] В частности, если $H_0 = H_0(x^*, \lambda^*)$, то $H_0 = \const$ на $t \in (0, T)$.
Детализация условий ПМП для двух типичных случаев: \[ U = \R, \qquad U = \R_+. \]
  1. 1 и 2 условия ПМП можно переписать как:
    1. $H_\nu \equiv 0, \qquad t \in [0, T]$
    2. $\dot H_\nu = 0$ при $t \in (\theta_i^*, \theta_{i+1}^*)$,
    где $\dot H_\nu = \dfrac{d H_\nu}{d t}$. Условия \[ \tag{12} \begin{aligned} H_\nu &\equiv 0, \qquad t \in [0, T] \\ \dot H_\nu &= 0, \qquad t \in (\theta_i^*, \theta_{i+1}^*) \end{aligned} \] являются условиями стационарности.
  2. Из 1 и 2 условия ПМП следует, что \[ H_\nu = \begin{cases} \leqslant 0, & t \in [0, T], \\ 0, & t \in S, \end{cases} \] где $S$ — все интервалы, где управление ненулевое ($u^*(t) \neq 0$): \[ (\theta_i^*, i = \overline{1, k}) \in S. \]
В формулировке ПМП предполагаем все функции гладкими по $x$. Условия 3 и 4 требуют дополнительной гладкости по времени. В противном случае условие 4 необходимо опустить, а в 3 условии дифференцируемость заменить на непрерывность.
ПМП распространяется на следующий случай. Пусть уравнения движения даны в виде \[ \dot x = f(t, x, v) + g(t, x) \nu, \] где $v = (v_1, \dots, v_r)$ — обычное кусочно-непрерывное управление, и \[ F_0 = F_0(t, x, v, \nu) = f_0(t, x, v) + g_0(t, x) \nu. \] Тогда условия 1-4 ПМП дополняются условием максимальности по $v$: \[ H_0(t, x^*, \lambda^*, v^*) = \max_{v \in V} H_0(t, x^*, \lambda^*, v), \] где $H_0 = \lambda f(x, t, v) - f_0(x, t, v)$ и \[ v \in V \subset \R^r, \] а $v^*(t)$ — оптимальное управление.
ПМП в общем является необходимым условием, но если задача является линейно выпуклой, то ПМП является и достаточным условием оптимальности (то есть если уравнения движения линейны по $x$, а функционал является выпуклым).
Ограничения на траекторию:
  1. Если граничные условия заданы в виде равенств, то условия трансверсальности не работают. То есть если \[ x(T) = x_\text{к}, \] то $\lambda(T)$ выбираются из выполнения этих граничных условий.
  2. Если ограничения имеют вид векторного неравенства: \[ x(T) \geqslant x_\text{к} \quad \text{либо} \quad x(T) \leqslant x_\text{к}, \] то появляются дополняющие условия нежёсткости: \[ \lambda(T) \gt 0 \quad \text{либо} \quad \lambda(T) \lt 0, \] причём \[ \dp{\lambda(T)}{(x^*(T) - x_\text{к})} = 0. \] Если условие активно ($x^*(T) = x_\text{к}$), то \[ \lambda(T) \lessgtr 0. \] Если же $x^*(T) \lessgtr x_\text{к}$, то \[ \lambda(T) = 0. \]
Посмотри примеры в тимсе.

§3. Примеры прикладных задач

1. Оптимизация расходов на рекламу в модели Эрроу-Нерлофа

Необходимо найти управление $u(t) \geqslant 0$, минимизирующее \[ J = \int\limits_{0}^{T} (u(t) - A \sqrt{x(t)}) dt \] для динамической системы \[ \tag{1} \dot x = u - bx, \qquad x(0) = x_0 \gt 0. \] Здесь

Тогда критерий качества $J$ имеет смысл суммарных издержек на рекламу. Минимизация $J$ равносильна максимизации прибыли от рекламы за период $[0, T]$.

Непрерывное выделение средств на рекламу моделируется поточечным ограничением $u(t) \leqslant c$. В этом случае задача классическая.

Рассмотрим единовременное выделение некоторой суммы $c$ на рекламу, распределение которой во времени оптимизируется.

В этом случае $U = \mathbb{R}_+$, и существование оптимального кусочно-непрерывного управления не гарантировано. Следовательно, сразу перейдём к импульсному управлению $\nu$.

Выпишем условия ПМП для рассматриваемого случая (см. 1 замечание к ПМП): \[ \tag{2} \begin{aligned} H &= \lambda (\nu - b x) - \nu + A \sqrt{x(t)}, \\ \dot \lambda &= \lambda b - \frac{A}{2 \sqrt{x(t)}}. \end{aligned} \] Тогда \[ H_\nu = \lambda - 1 = \begin{cases} \leqslant 0, & t \in [0, T] \setminus S, \\ = 0 & t \in S. \end{cases} \] Так как $\nu$ не входит в сопряжённое уравнение (уравнение Эйлера-Лагранжа), то $\lambda(t)$ — кусочно-гладкая функция.

Проанализируем необходимые условия в обратном времени, то есть от $T$ к $t = 0$, используя граничные условия для $\lambda$.

Так как $\varphi_0 \equiv 0$, то \[ \varphi_0(T) = 0, \implies \lambda(T) = 0, \implies H_\nu(T) = \lambda(T) - 1 = -1 \lt 0, \] поэтому, в силу условия непрерывности $H_\nu$ на $[0, T]$ (условие 3 ПМП), $H_\nu(t) \lt 0$ по крайней мере в некоторой левой полуокрестности $O(T)$ точки $T$. Следовательно, $\nu^* = 0$ в $O(T)$.

Составим производную функции переключения по $t$: \[ \tag{3} \dot H_\nu = \dot \lambda = b \lambda^* - \frac{A}{2 \sqrt{x^*(t)}}. \] Так как $\lambda^*(T) = 0$, то \[ \dot H_\nu \approx - \frac{A}{2 \sqrt{x^*(t)}} \lt 0 \] в окрестности $O(T)$. Следовательно, функция переключения строго убывает по $t$ в полуокрестности $O(T)$ вдоль оптимального процесса. Это значит, что левее точки $T$ может оказаться момент времени $t_1 \in (0, T)$, где $H_\nu = 0$, то есть возможный момент импульса или конец интервала, где $u^*(t) \gt 0$.

Если $T$ достаточно мало, то такого момента нет, и оптимальное управление $\nu^* = 0$ на $[0, T]$. Это соответствует случаю, когда вложения в рекламу вообще не выгодны, поскольку их эффективность не успевает проявиться.

Пусть $T$ теперь достаточно велико, и момент $t_1 \in (0, T)$ существует. Тогда \[ H_\nu(t_1) = \lambda(t_1) - 1 = 0, \implies \lambda(t_1) = 1. \]

Покажем, что $t_1$ не может быть моментом импульса. Действительно, если \[ H_\nu(t_1) = 0, \] то $H_\nu$ имеет в $t_1$ локальный максимум (так как $H_\nu \lt 0$ в $O(T)$). Следовательно, верны неравенства: \[ \dot H_\nu(t_{1-}) \geqslant 0, \qquad \dot H_\nu(t_{1+}) \leqslant 0. \] Из (3) и непрерывности $\lambda(t)$ получаем эквивалентное неравенство: \[ \frac{A}{2 \sqrt{x(t_{1-})}} \lt \frac{A}{2 \sqrt{x(t_{1+})}}, \] откуда \[ x(t_{1-}) \gt x(t_{1+}). \] Но из условия скачка для системы (1) получаем \[ \dv{z}{\tau} = c, \implies z(\tau) = c \tau + x(t_{1-}), \] но $z(1) = x(t_{1+})$, поэтому \[ x(t_{1+}) = c + x(t_{1-}), \implies x(t_{1+}) - x(t_{1-}) = c \gt 0 \implies x(t_{1+}) \gt x(t_{1-}). \] Значит, момент времени $t_1$ не удовлетворяет условию скачка, поэтому он является правым концом интервала $\Delta$, на котором $u^*(t) \gt 0$. Тогда из первого замечания к ПМП следует, что \[ H_\nu \equiv 0 \quad \text{и} \quad \dot H_\nu = 0 \quad \text{на} \quad \Delta, \] и из (3) получаем \[ \tag{4} \lambda^* = 1, \implies x^* = \frac{A^2}{4 b^2} \quad \text{на} \quad \Delta. \] Подстановка этого значения $x^*$ в уравнение (1) позволяет найти управление на особом участке, где $u^*(t) \gt 0$.

Продолжая рассуждения в обратном времени, обозначим левый конец интервала $\Delta$ как $t_0$. Поскольку значение $x^*$ из (4) может совпасть с $x_0$ только случайно, то попаданию на особый участок $\Delta$ должно предшествовать управление, отличное от особого.

Опять существуют две возможности: $t_0$ — момент импульса и $t_0$ — момент непрерывного попадания траектории на особый участок.

Первый вариант реализуется, если \[ x_0 \lt \frac{A^2}{4 b^2}, \qquad t_0 = 0. \] Тогда величина импульса \[ c^* = \frac{A^2}{4 b^2} - x_0 \] находится из условия скачка.

Второй вариант возможен если \[ x_0 \gt \frac{A^2}{4 b^2}. \] Это означает высокий начальный уровень расположенности покупателей, поддержание которого экономически не выгодно. В этом случае моменту $t_0$ предшествует управление $u^* = 0$ и экспоненциальное уменьшение $x^*(t)$: \[ \dot x = -b \implies x^* = e^{-bt}. \]

2. Оптимизация рекламных расходов в модели Видала-Вулфа
2.1. Постановка задачи

Авторы модели исходят из допущения, что изменение объёма продаж обусловлено двумя факторами:

  1. часть информированного рекламой населения либо забывает об этом товаре, либо отказывается от услуг фирмы, либо покидает рынок с постоянным темпом, если фирма не вкладывается в рекламу;
  2. увеличение темпа продаж пропорционально текущим вложениям в рекламу $u(t)$, а также части рынка, пока не охваченного, но потенциального.

Логика математического моделирования процесса такова.

Таким образом, составленная мат. модель имеет вид \[ \tag{5} \dot s(t) = k u(t) \paren{ 1 - \frac{s(t)}{M} } - b s(t), \qquad s(0) = s_0, \] где $k = \const$ — коэффициент пропорциональности.

Заметим, что при достаточно большом уровне насыщения $M \to \infty$ уравнение (5) принимает вид \[ \dot s(t) = k u(t) - b s(t), \qquad s(0) = s_0, \] то есть структуру линейного уравнения модели Эрроу-Нерлофа.

Итак, рассмотрим (5). Сделаем замену переменной: \[ x = \frac{s(t)}{M}, \] где $x$ — доля фактического объёма продаж от потенциально возможного. Положим \[ a = \frac{k}{M}, \] тогда (5) примет вид \[ \dot x(t) = a u(t) \paren{ 1 - x(t) } - b x(t), \qquad x(0) = x_0 \in (0, 1). \] Для оптимизации рекламных расходов максимизируем дисконтированную прибыль за период времени $T$, то есть будем рассматривать задачу Лагранжа \[ \tag{6} J = \int\limits_{0}^{T} e^{-rt} \paren{ p x(t) - u(t) } dt \to \sup \] при ограничениях \[ \tag{7} u(t) \geqslant 0, \qquad \int\limits_{0}^{T} u(t) dt \leqslant R. \] Здесь

2.2. Исследование модели

Приведём поставленную задачу (6) к задаче Майера. Для этого введём дополнительные фазовые переменные: \[ y(t) = \int\limits_{0}^{t} e^{-r \tau} \paren{ p x(\tau) - u(\tau) } d\tau, \qquad \gamma(t) = \int\limits_{0}^{t} u(\tau) d\tau. \] Тогда исходная задача примет вид \[ \tag{8} \begin{aligned} \dot x &= -b x(t) + a u(t) (1 - x(t)), & x(0) &= x_0 \in (0, 1) \\ \dot y &= e^{-r t} (p x(t) - u(t)), & y(0) &= 0, \\ \dot \gamma &= u(t), & \gamma(0) &= 0, \\ && \gamma(T) &\leqslant R. \end{aligned} \] Целевой функционал: \[ J(u) = -y(T) \to \inf, \] огранчение на управление: $u(t) \geqslant 0$.

Задачу (8) будем рассматривать в импульсной постановке, то есть $u(t) \to \nu(t)$.

Составим гамильтониан: \[ H = \lambda_1 \left[ a \nu(t) (1 - x(t)) - b x(t) \right] + \lambda_2 e^{-rt} \paren{ p x(t) - \nu(t) } + \lambda_3 \nu(t). \] Согласно ПМП $\lambda \neq 0$. Условия трансверсальности, согласно (8), будут работать только для $\lambda_1(T), \lambda_2(T)$. Значение же $\lambda_3(T)$ должно удовлетворять условию дополняющей нежёсткости: \[ (9) \lambda_3(T) \leqslant 0, \qquad (\lambda_3(T), \gamma(T) - R) = 0, \] то есть если ограничение на $\gamma(T)$ активно и $\gamma(T) = R$, то $\lambda_3(T) \lt 0$; если же $\gamma(T) \lt R$, то $\lambda_3(T) = 0$.

Пусть $\beta = \const \geqslant 0$, тогда $\lambda_3(T) = - \beta$. В этом случае уравнения Эйлера-Лагранжа следующие: \[ \begin{aligned} \dot \lambda_1 &= - \pd{H}{x} = \lambda_1(a \nu(t) + b) - \lambda_2 e^{-rt} p, & \lambda_1(T) &= - \pd{\varphi_0(T)}{x} = 0, \\ \dot \lambda_2 &= - \pd{H}{y} = 0, & \lambda_2(T) &= - \pd{\varphi_0(T)}{y} = 1, \implies \lambda_2 \equiv 1, \\ \dot \lambda_2 &= - \pd{H}{\gamma} = 0, & \lambda_3(T) &= -\beta. \end{aligned} \] Ограничение на $\nu$ такое: $\nu \in U = \mathbb{R}_+$, поэтому из ПМП следует, что (см. замечание 1) \[ \begin{aligned} H_\nu &= \lambda_1 a(1 - x) - \lambda_2 e^{-rt} + \lambda_3 = \\ &= \lambda_1 a(1 - x) - e^{-rt} - \beta = \\ &= \begin{cases} \leqslant 0, & \forall t, \\ = 0 & t \in S(\nu). \end{cases} \end{aligned} \]

Исследуем структуру оптимального управления.

  1. Рассмотрим $t = T$, тогда \[ H_\nu = e^{-rt} - \beta \lt 0, \] следовательно, нет конечного импульса и $H_\nu \lt 0$ слева от $T$ в силу непрерывности $H_\nu$. Это означает, что оптимальное управление равно нулю в некоторой левой полуокрестности $O(T)$ точки $T$.
  2. Выведем условие невыгодности инвестиций. Пусть $\nu(t) \equiv 0$, тогда \[ \dot \gamma = 0 \implies \gamma = \const, \] и из условия $\gamma(0) = 0$ следует, что \[ \gamma(t) \equiv 0 \lt R. \] Значит, ограничение на $\gamma(T)$ не активно, поэтому $\lambda_3(T) = 0$, но из системы (10) следует, что в этом случае $\lambda_3(T) \equiv 0$.

    Тогда из системы (10) следует, что \[ \dot \lambda_1 = \lambda_1 b - e^{-rt} p, \] значит, \[ \lambda_1(t) = C_1 e^{bt} + \frac{p}{b + r} e^{-rt}. \] Из условия $\lambda_1(T) = 0$ следует, что \[ C_1 e^{bT} + \frac{p}{b + r} e^{-rT} = 0, \implies C_1 = - \frac{p}{b + r} e^{-rT - bT}. \] Значит, \[ \lambda_1(t) =\frac{p}{b + r} \left[ e^{-rt} - e^{bt - (b + r) T} \right]. \] Так как в случае $\nu(t) \equiv 0$ из (8) следует, что \[ x(t) = x_0 e^{-bt}, \] и, учитывая, что $\lambda_3(t) \equiv 0 \implies \beta = 0$, из условия \[ H_\nu \lt 0 \] следует, что \[ a \lambda_1(t) (1 - x_0 e^{-bt}) \lt e^{-rt}, \qquad t \in [0, T]. \]

  3. Пусть $\nu(t) \not\equiv 0$. Рассмотрим особый интервал $\Delta$, на котором обычная компонента управления $u(t) \gt 0$. Тогда (см. замечание 1 к ПМП) \[ H_\nu = 0 \quad \text{и} \quad \dot H_\nu = 0 \quad \text{на} \quad \Delta, \] то есть \[ H_\nu = 0, \implies \lambda_1 a(1 - x) - e^{-rt} - \beta = 0, \] значит, \[ \lambda_1 = \frac{e^{-rt} + \beta}{a(1 - x)}. \] Тогда из условия $\dot H_\nu = 0$ получаем квадратное уравнение относительно $x(t)$, откуда \[ \tag{11} x^* = \frac{2 B + D - \sqrt{D^2 + 4 AB}}{2 B}, \] где \[ A = b\paren{ e^{-rt} + \beta }, \qquad B = ape^{-rt}, \qquad D = re^{-rt}. \] Значение $x^*$ допустимо (то есть $x^* \in (0,1)$) при выполнении условия \[ b + r \lt a p. \] Магистральное управление $u^*$ находится из уравнения \[ u^*(t) = \frac{b x^*}{a(1 - x^*)}. \]
  4. Непрерывное (при $x_0 \gt x^*$) и импульсное (при $x_0 \lt x^*$) попадание на магистраль определяется как в предыдущей задаче. Величина $c^*$ начального импульса ищется из предельной системы исходного уравнения \[ \dv{z}{\tau} = a c^* (1 - z), \qquad z(0) = x_0, \qquad z(1) = x^* \] и равна \[ c^* = \frac{1}{a} \ln \frac{1 - x_0}{1 - x^*}. \]
  5. Магистраль (11) зависит от $\beta \geqslant 0$. Чтобы конкретизировать полученное решение, предположим, что ресурсы не дефицитны, то есть: Если полученное таким путём решение не удовлетворяет всем условиям допустимости и ПМП, то это означает полное расходование всех средств, следовательно, $\gamma(T) = R$, и необходимо подбирать $\beta \gt 0$ за счёт варьирования времени пребывания на магистрали.