<jupyter_start><jupyter_text>Laboratório 9: Populações de Urso Referente ao capítulo 15 O cenário neste laboratório é motivado pelo *Great Smoky Mountain National Park*. É considerad a população de ursos em um parque genérico em áreas populadas por humanos, tal que essas regiões apresentem bordas compartilhadas. A emigração e imigração de ursos é baseada na conectividade entre caminhos. |Parâmetro|Significado||---------|-----------||$r$|taxa de crecimento de ursos||$K$|capacidade de carga||$m_p$|proporção da fronteira do parque conectada com a floresta||$m_f$|proporção da fronteira da floresta conectada com o parque|Seja $P(t), F(t)$ e $O(t)$ a densidade de ursos no parque, na floresta e nas regiões de fora. Os controles serão $u_p$ e $u_f$, que são as taxas de caça de ursos no parque e na floresta, respectivamente. Atualmente a caça em parques não é legalizada, mas o modelo considera que ocorrem também, além da caça na floresta. Queremos minimizar o número de ursos nas áreas populadas por humanos. Assim, o problema pode ser descrito pelo seguinte modelo. $$ \min_{u_p, u_f} \int_0^T O(t) + c_pu_p(t)^2 + c_fu_f(t)^2 dt $$$$\text{sujeito a  }P'(t) = rP(t) - \frac{r}{K}P(t)^2 + \frac{m_fr}{K}\left(1 - \frac{P(t)}{K}\right)F(t)^2 - u_p(t)P(t),$$$$F'(t) = rF(t) - \frac{r}{K}F(t)^2 + \frac{m_pr}{K}\left(1 - \frac{F(t)}{K}\right)P(t)^2 - u_f(t)F(t),$$$$O'(t) = r(1 - m_p)\frac{P(t)^2}{K} + r(1 - m_f)\frac{F(t)^2}{K} + \frac{m_fr}{K^2}P(t)F(t)^2 + \frac{m_pr}{K^2}P(t)^2F(t),$$$$P(0) = P_0 \geq 0, F(0) = F_0 \geq 0, O(0) = O_0 \geq 0, $$$$0 \leq u_p(t) \leq 1, 0 \leq u_f(t) \leq 1 $$ O termo $\frac{r}{K}P(t)^2$ representa a emigração densidade-depentente. O terceiro termo na primeira equação é a imigração da floresta densidade-dependente e $u_pP$ o nível de caça. O primeiro termo da terceira equação é a emigração do parque para regiões fora do parque, enquanto o segundo é a emigração da floresta para outras regiões.  Importanto as bibliotecas<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import sympy as sp

import sys  
sys.path.insert(0, '../pyscripts/')

from optimal_control_class import OptimalControl<jupyter_output><empty_output><jupyter_text>Condições NecessáriasVamos utilizar a bilbioteca `sympy` para nos ajudar a escrever as condições necessárias.<jupyter_code>P, F, O, uf, up, l1, l2, l3, cp, cf, r, K, mf, mp, t = sp.symbols('P F O uf up l1 l2 l3 cp cf r K mf mp t')<jupyter_output><empty_output><jupyter_text>O Hamiltoniano\begin{multline}H(t, \vec{x}, \vec{u}, \vec{\lambda}) = O(t) + c_pu_p(t)^2 + c_fu_f(t)^2 +  \\ \lambda_1(t)\left[rP(t) - \frac{r}{K}P(t)^2 + \frac{m_fr}{K}\left(1 - \frac{P(t)}{K}\right)F(t)^2 - u_p(t)P(t)\right] + \\ \lambda_2(t)\left[rF(t) - \frac{r}{K}F(t)^2 + \frac{m_pr}{K}\left(1 - \frac{F(t)}{K}\right)P(t)^2 - u_f(t)F(t)\right] + \\ \lambda_3(t)\left[r(1 - m_p)\frac{P(t)^2}{K} + r(1 - m_f)\frac{F(t)^2}{K} + \frac{m_fr}{K^2}P(t)F(t)^2 + \frac{m_pr}{K^2}P(t)^2F(t)\right]\end{multline}<jupyter_code>H =  O + cp*up**2 + cf*uf**2 + l1*(r*P - r*P**2/K + mf*(r/K)*(1 - P/K)*F**2 - up*P) + \
                               l2*(r*F - r*F**2/K + mp*(r/K)*(1 - F/K)*P**2 - uf*F) + \
                               l3*(r*(1 - mp)*P**2/K + r*(1 - mf)*F**2/K + mf*r*P*F**2/K**2 + mp*r*P**2*F/K**2)<jupyter_output><empty_output><jupyter_text>Condições do Estado<jupyter_code>print(sp.diff(H,l1), '\n')
print(sp.diff(H,l2), '\n')
print(sp.diff(H,l3))<jupyter_output>F**2*mf*r*(1 - P/K)/K + P*r - P*up - P**2*r/K 

-F**2*r/K + F*r - F*uf + P**2*mp*r*(-F/K + 1)/K 

F**2*r*(1 - mf)/K + F**2*P*mf*r/K**2 + F*P**2*mp*r/K**2 + P**2*r*(1 - mp)/K<jupyter_text>Equações Adjuntas<jupyter_code>dl1 = (-1)*sp.diff(H,P)
dl2 = (-1)*sp.diff(H,F)
dl3 = (-1)*sp.diff(H,O)
print(dl1, '\n')
print(dl2, '\n')
print(dl3)<jupyter_output>-l1*(-F**2*mf*r/K**2 + r - up - 2*P*r/K) - l3*(F**2*mf*r/K**2 + 2*F*P*mp*r/K**2 + 2*P*r*(1 - mp)/K) - 2*P*l2*mp*r*(-F/K + 1)/K 

-2*F*l1*mf*r*(1 - P/K)/K - l2*(-2*F*r/K + r - uf - P**2*mp*r/K**2) - l3*(2*F*r*(1 - mf)/K + 2*F*P*mf*r/K**2 + P**2*mp*r/K**2) 

-1<jupyter_text>Condições de Otimalidade<jupyter_code>du1 = sp.diff(H,up)
du2 = sp.diff(H,uf)
print(du1, '\n')
print(du2)<jupyter_output>-P*l1 + 2*cp*up 

-F*l2 + 2*cf*uf<jupyter_text>$$ \frac{\partial H}{\partial u_p} > 0 \implies u_p^*(t) = 0 \implies \frac{P\lambda_1}{2c_p} < 0$$$$ \frac{\partial H}{\partial u_p}  1$$$$ \frac{\partial H}{\partial u_p} = 0 \implies 0 \leq u_p^*(t) = \frac{P\lambda_1}{2c_p} \leq 1$$$$u_p^* = \max(0, \min(1, \frac{P\lambda_1}{2c_p}))$$$$ \frac{\partial H}{\partial u_f} > 0 \implies u_f^*(t) = 0 \implies \frac{F\lambda_2}{2c_f} < 0$$$$ \frac{\partial H}{\partial u_f}  1$$$$ \frac{\partial H}{\partial u_f} = 0 \implies 0 \leq u_f^*(t) = \frac{F\lambda_2}{2c_f} \leq 1$$$$u_f^* = \max(0, \min(1, \frac{F\lambda_2}{2c_f}))$$  Aplicando a classe ao exemplo Nesse laboratório, o leitor é convidado a escrever as equações e realizar os testes por conta. Você deve mudar a célula de Raw para Code e escrever as funções como fizemos até então.<jupyter_code>parameters = {'r': None, 'K': None, 'mp': None, 'mf': None, 'cp': None, 'cf': None}

# Lembre: P: x[0], F: x[1], O: x[2], up: u[0], uf: u[1]
diff_state = lambda t, x, u, par: np.array([
    x[1]**2*par['mf']*par['r']*(1 - x[0]/par['K'])/par['K'] + x[0]*par['r'] - \
        x[0]*u[0] - x[0]**2*par['r']/par['K'], 
    -x[1]**2*par['r']/par['K'] + x[1]*par['r'] - x[1]*u[1] + \
        x[0]**2*par['mp']*par['r']*(-x[1]/par['K'] + 1)/par['K'], 
    x[1]**2*par['r']*(1 - par['mf'])/par['K'] + x[1]**2*x[0]*par['mf']*par['r']/par['K']**2 + \
        x[1]*x[0]**2*par['mp']*par['r']/par['K']**2 + x[0]**2*par['r']*(1 - par['mp'])/par['K']
])

diff_lambda = lambda t, x, u, l, par: np.array([
    -l[0]*(-x[1]**2*par['mf']*par['r']/par['K']**2 + par['r'] - u[0] - 2*x[0]*par['r']/par['K']) - \
        l[2]*(x[1]**2*par['mf']*par['r']/par['K']**2 + 2*x[1]*x[0]*par['mp']*par['r']/par['K']**2 + \
        2*x[0]*par['r']*(1 - par['mp'])/par['K'])-2*x[0]*l[1]*par['mp']*par['r']*(1-x[1]/par['K'])/par['K'], 
    -2*x[1]*l[0]*par['mf']*par['r']*(1 - x[0]/par['K'])/par['K'] - \
        l[1]*(-2*x[1]*par['r']/par['K'] + par['r'] - u[1] - x[0]**2*par['mp']*par['r']/par['K']**2) - \
        l[2]*(2*x[1]*par['r']*(1 - par['mf'])/par['K'] + 2*x[0]*x[1]*par['mf']*par['r']/par['K']**2 + \
        x[0]**2*par['mp']*par['r']/par['K']**2), 
    -1
])

update_u = lambda t, x, lambda_, par: np.array([
    np.maximum(0.0, np.minimum(1.0, 0.5*x[0]*lambda_[0]/par['cp'])),
    np.maximum(0.0, np.minimum(1.0, 0.5*x[1]*lambda_[1]/par['cf']))
])
problem = OptimalControl(diff_state, diff_lambda, update_u, 
                         n_controls = 2, n_states = 3, 
                         bounds = [(0,1), (0,1)])
x0 = np.array([0.4, 0.2, 0])
T = 25
parameters['r'] = 0.1
parameters['K'] = 0.75
parameters['mp'] = 0.5
parameters['mf'] = 0.5
parameters['cp'] = 10000
parameters['cf'] = 10
t,x,u,lambda_ = problem.solve(x0, T, parameters)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.26519572953736653, 'fraction_numerical': 0.04313167259786477, 'mean_word_length': 4.072924187725632, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '222704', 'license': 'MIT', 'n_tokens_mistral': 3505, 'n_tokens_neox': 3410, 'n_words': 640, 'path': 'notebooks/.ipynb_checkpoints/Laboratory9-checkpoint.ipynb', 'repo_name': 'lucasmoschen/Optimal_Control_Biological', 'size': '7025'}