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

[re3burn]

Добрый день!
При разборе статьи по оптимальному управлению моделью иммунологии ВИЧ возникло пару вопросов. Оптимальное управление строится с использованием принципа максимума.
Модель динамики:

$$ \begin{cases} \dot{S_{1}}=\mu_{1}-d_{1}S_{1}-(1-u_{1})k_{1}VS_{1}\\ \dot{S_{2}}=\mu_{2}-d_{2}S_{2}-(1-fu_{1})k_{2}VS_{2}\\ \dot{I_{1}}=(1-u_{1})k_{1}VS_{1}-\delta I_{1}-m_{1}EI_{1}\\ \dot{I_{2}}=(1-fu_{1})k_{2}VS_{2}-\delta I_{2}-m_{2}EI_{2}\\ \dot{V}=(1-u_{2})N\delta (I_{1}+I_{2})-cV-[(1-u_{1})\rho_{1} k_{1}S_{1}+(1-fu_{1})\rho_{2} k_{2}S_{2}]V\\ \dot{E}=\mu_{1}+\frac {b_{E}(I_{1}+I_{2})} {(I_{1}+I_{2})K_{b}} E - \frac {d_{E}(I_{1}+I_{2})} {(I_{1}+I_{2})K_{d}} E - \delta_{E} E\\ \end{cases} $$

Минимизируемый функционал качества:

$$ \mathcal{J}(u_1,u_2,t)=\frac {1} {2} \left \lbrace \int\limits_0^T (R_1u_1^2+R_2u_2^2)dt + Q(V-\tilde{V})^2+S(E-\tilde{E})^2+PT^2 \right \rbrace $$

От времени t зависят переменные состояния $ \bar{x}(t)=[S_{1}(t), S_{2}(t), I_{1}(t), I_{2}(t), V(t), E(t)]$ и переменные управления $ \bar{u}(t)=[u_{1}(t), u_{2}(t)] $. Остальные переменные - константы.
Левый конец закреплен $\bar{x}(0)=\bar{x}_0$, правый - свободный. На управления наложены ограничения $0\leq a_1\leq u_1\leq b_1\leq 1, 0\leq a_2\leq u_2\leq b_2\leq 1$
Итого, в стандартных обозначениях задача оптимального управления принимает вид:

$$ \begin{cases} \frac{d\bar{x}}{dt} = \bar{f}(\bar{x}(t),\bar{u}(t),t) \\ \bar{u}(t)\in U\\ \bar{x}(0) = \bar{x}_0\\ \mathcal{J}(\bar{u},t)=\int\limits_0^T f_{0}(\bar{x}(t),\bar{u}(t),t)dt + F(\bar{x}(T),T) \to \min\\ \end{cases} $$

Гамильтониан системы:

$$ \mathcal{H}(\bar{x}(t),\bar{u}(t),\bar{\psi}(t),t,\lambda_{0})=-\lambda_0 f_0(\bar{x}(t),\bar{u}(t),t)+\sum\limits_{j=1}^{6} \psi_j(t)f_j(\bar{x}(t),\bar{u}(t),t)=$$ $$=[\text{для задачи со свободным правым концом } \lambda_{0}=1]=$$ $$=-\frac {1} {2} (R_1u_1^2+R_2u_2^2)+\psi_1(t)\left[\mu_{1}-d_{1}S_{1}-(1-u_{1})k_{1}VS_{1} \right]+$$ $$+\psi_2(t)\left[\mu_{2}-d_{2}S_{2}-(1-fu_{1})k_{2}VS_{2} \right]+$$ $$+\psi_3(t)\left[(1-u_{1})k_{1}VS_{1}-\delta I_{1}-m_{1}EI_{1} \right]+$$ $$+\psi_4(t)\left[(1-fu_{1})k_{2}VS_{2}-\delta I_{2}-m_{2}EI_{2} \right]+$$ $$+\psi_5(t)\left[(1-u_{2})N\delta (I_{1}+I_{2})-cV-[(1-u_{1})\rho_{1} k_{1}S_{1}+(1-fu_{1})\rho_{2} k_{2}S_{2}]V \right]+$$ $$+\psi_6\left[\mu_{1}+\frac {b_{E}(I_{1}+I_{2})} {(I_{1}+I_{2})K_{b}} E - \frac {d_{E}(I_{1}+I_{2})} {(I_{1}+I_{2})K_{d}} E - \delta_{E} E\right] $$

Для решаемой задачи сопряженная система:

$$ \begin{cases} \dot{\psi_1}=-\frac {\partial \mathcal{H}}{\partial S_1}=-\left \lbrace \psi_1[-d_1-(1-u_1)k_1V]+\psi_3(1-u_1)k_1V-\psi_5(1-u_1)\rho_1k_1V \right \rbrace\\ \dot{\psi_2}=-\frac {\partial \mathcal{H}}{\partial S_2}=-\left \lbrace \psi_2[-d_2-(1-fu_1)k_2V]+\psi_4(1-fu_1)k_2V-\psi_5(1-fu_1)\rho_2k_2V \right \rbrace\\ \dot{\psi_3}=-\frac {\partial \mathcal{H}}{\partial I_1}=-\left \lbrace \psi_3(-\delta-m_1E)+\psi_5(1-u_2)N\delta+\psi_6 \left(\frac {b_EEK_b} {(I_1+I_2+K_b)^2} - \frac {d_EEK_d} {(I_1+I_2+K_d)^2} \right) \right \rbrace\\ \dot{\psi_4}=-\frac {\partial \mathcal{H}}{\partial I_2}=-\left \lbrace \psi_3(-\delta-m_2E)+\psi_5(1-u_2)N\delta+\psi_6 \left(\frac {b_EEK_b} {(I_1+I_2+K_b)^2} - \frac {d_EEK_d} {(I_1+I_2+K_d)^2} \right) \right \rbrace\\ \dot{\psi_5}=-\frac {\partial \mathcal{H}}{\partial V}=-\lbrace \-\psi_1(1-u_1)k_1S_1-\psi_2(1-fu_1)k_2S_2+\psi_3(1-u_1)k_1S_1+\psi_4(1-fu_1)k_2S_2-\\ \quad\quad\quad\quad\quad-\psi_5[c+(1-u_1)\rho_1k_1S_1+(1-fu_1)\rho_2k_2S_2] \rbrace\\ \dot{\psi_6}=-\frac {\partial H}{\partial E}=-\left \lbrace -\psi_3m_1I_1-\psi_3m_2I_2+\psi_6\left[\frac {b_{E}(I_{1}+I_{2})} {(I_{1}+I_{2})K_{b}} - \frac {d_{E}(I_{1}+I_{2})} {(I_{1}+I_{2})K_{d}} - \delta_{E} \right] \right \rbrace\\ \end{cases} $$

Далее авторы статьи строят расширенный штрафными функциями гамильтониан(лагранжиан):

$$ \mathcal{L}(\bar{x}(t),\bar{u}(t),\bar{\psi}(t),t)=\mathcal{H}(\bar{x}(t),\bar{u}(t),\bar{\psi}(t),t)-$$ $$\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad -\omega_{11}(t)(u_1(t)-a_1)-\omega_{12}(t)(b_1-u_1(t))-\omega_{21}(t)(u_2(t)-a_2)-\omega_{22}(t)(b_2-u_2(t)) $$

Штрафные слагаемые $\omega_{ij}(t)\geq 0$ и

$$ \omega_{11}(t)(u_1(t)-a_1)=\omega_{12}(t)(b_1-u_1(t))=0 \text{,если } u_1(t)=u_1^{*}$$ $$\omega_{21}(t)(u_2(t)-a_2)=\omega_{22}(t)(b_2-u_2(t))=0 \text{,если } u_2(t)=u_2^{*} $$

Далее приводятся альтернативные выражения для сопряженной системы $\dot{\psi}_j(t)$:

$$ \dot{\psi}_1=\frac {\partial{\mathcal{L}}}{\partial{S_1}},\quad\quad \dot{\psi}_2=\frac {\partial{\mathcal{L}}}{\partial{S_2}},\quad\quad \dot{\psi}_3=\frac {\partial{\mathcal{L}}}{\partial{I_1}},\quad\quad \dot{\psi}_4=\frac {\partial{\mathcal{L}}}{\partial{I_2}},\quad\quad \dot{\psi}_5=\frac {\partial{\mathcal{L}}}{\partial{V}},\quad\quad \dot{\psi}_6=\frac {\partial{\mathcal{L}}}{\partial{E}}. $$

Условия трансверсальности:
1) в конечный момент времени
$t=T:\text{ }\psi_j(T)=\frac{\partial{F(x_j(T),T)}}{\partial{x_j}}\Rightarrow \bar{\psi}(T^{*})=\left(0, 0, 0, 0, Q(V(T^{*})-\tilde{V}), S(E(T^{*})-\tilde{E}) \right)$
2) из условия равенства первой вариации функционала нулю:

$$ \lim_{\delta\rightarrow 0}{\frac{\mathcal{J}(u_1^*,u_2^*,T^*+\delta)-\mathcal{J}(u_1^*,u_2^*,T^*)}{\delta}}\equiv 0 \Leftrightarrow$$ $$\Leftrightarrow \frac 1 2 \left ( R_1u_1^*^2 + R_2u_2^*^2\right) +\psi_5(T^*)\dot{V}^*(T^*)+\psi_6(T^*)\dot{E}^*(T^*)+PT^* $$

Выражения для оптимальных управлений $u_1^*$ и $u_2^*$:

$$ u_1^*=\max\left (a_1, \min \left( b_1, -\frac {-(\psi_1-\psi_3+\rho_1\psi_5)k_1VS_1 -(\psi_2-\psi_4+\rho_2\psi_5)fk_2VS_2} {R_1}\right)\right)$$ $$u_2^*=\max\left (a_2, \min \left( b_2, -\frac {\psi_5N\delta(I_1+I_2)} {R_2}\right)\right) $$

У меня возникли следующие вопросы:
1) у авторов статьи в выражениях для гамильтониана коэффициент 1/2 входит со знаком "+" и знаки в выражениях для оптимальных управлений противоположны (правильно ли у меня выбрано $\lambda_0=1$) ?
2) для чего нужны альтернативные выражения $\dot{\psi}_j=\frac {\partial{\mathcal{L}}}{\partial{x_j}}$ для сопряженной системы, если сопряженная система уже получена ранее?
3) необходимы ли условия трансверсальности 2) (равенство нулю первой вариации функционала) для решения задачи? Каким образом у авторов в этом уравнении появились $\psi_5(T^*)$ и $\psi_6(T^*)$?
4) как собственно получить численное решение задачи (промоделировать систему)?

С уважением, Алексей