PP

environment.yml

channels:
  - defaults
dependencies:
  - ipython
  - numpy
  - scipy
  - matplotlib
  - notebook
  - console_shortcut
  - ipywidgets
  - pip
  - pip:
     - gekko

Predator_Prey_Intro.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "01d7bb8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact, widgets, Layout\n",
    "from gekko import GEKKO\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9fe85e13",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "02c0924a3e8b41b8aab6a8c04da8d655",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.5, description='Initial Prey Population $\\\\quad$', layout=Layout(wid…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_pp(a, b)>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def model(x):\n",
    "    return np.array([x[0]-x[0]*x[1]-0.4*x[0],\n",
    "                     -x[1]+x[0]*x[1]-0.2*x[1]])\n",
    "def plot_pp(a,b):\n",
    "    t0 = 0\n",
    "    tmax = 12\n",
    "    sol = solve_ivp(lambda t, x:model(x), [t0, tmax], np.array([a,b]), t_eval=np.linspace(t0, tmax, 100))\n",
    "    plt.subplot(1,2,1)\n",
    "    plt.plot(sol.t,sol.y[0], label = 'Prey')\n",
    "    plt.plot(sol.t,sol.y[1], label = 'Predator')\n",
    "    plt.xlabel(\"time\")\n",
    "    plt.ylabel(\"population\")\n",
    "    plt.title(\"Solution of the ODE system\")\n",
    "    plt.legend()\n",
    "    plt.subplot(1,2,2)\n",
    "    plt.plot(sol.y[1],sol.y[0], label = 'Prey', color = \"r\")\n",
    "    plt.rcParams[\"figure.figsize\"]=10,5\n",
    "    plt.title(\"Phase portrait of the solution\")\n",
    "    plt.xlabel(\"predator population\")\n",
    "    plt.ylabel(\"prey population\")\n",
    "\n",
    "style = {'description_width': 'initial'}\n",
    "layout = Layout(width = \"400px\")\n",
    "slid1 = widgets.FloatSlider(description = 'Initial Prey Population $a$ $\\quad$',min=0, max=1, step=0.1, value=0.5,style = style, layout = layout)\n",
    "slid2 = widgets.FloatSlider(description = 'Initial Predator Population $b$',min=0, max=1, step=0.1, value=0.7, style = style, layout = layout)\n",
    "interact(plot_pp, a= slid1 ,b=slid2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f79e3f70",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "be27af2d939d44179af7e3684cf6916c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.5, description='Initial Prey Population $\\\\quad$', layout=Layout(wid…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_pp(a, b)>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def u(t):\n",
    "    if 0 <= t <= 2:\n",
    "        return 0\n",
    "    elif 2 < t <= 4.5:\n",
    "        return 1\n",
    "    elif 4.5 < t <= 7.5:\n",
    "        return 0.5\n",
    "    else:\n",
    "        return 0\n",
    "vu = np.vectorize(u, otypes=[np.float64])\n",
    "\n",
    "def model_control_ex(t, x):\n",
    "    return np.array([x[0] - x[0] * x[1] - 0.4 * x[0] * u(t), -x[1] + x[0] * x[1] - 0.2 * x[1] * u(t)])\n",
    "\n",
    "\n",
    "def plot_pp(a, b):\n",
    "    t0 = 0\n",
    "    tmax = 12\n",
    "    sol = solve_ivp(lambda t, x: model_control_ex(t, x), [t0, tmax], np.array([a, b]),\n",
    "                    t_eval=np.linspace(t0, tmax, 100))\n",
    "    plt.subplot(1, 2, 1)\n",
    "    plt.plot(sol.t, sol.y[0], label='Prey')\n",
    "    plt.plot(sol.t, sol.y[1], label='Predator')\n",
    "    plt.xlabel(\"time\")\n",
    "    plt.ylabel(\"population\")\n",
    "    plt.title(\"Solution of the ODE system\")\n",
    "    plt.legend()\n",
    "    plt.subplot(1, 2, 2)\n",
    "    plt.plot(sol.y[1], sol.y[0], label='Prey', color=\"r\")\n",
    "    plt.rcParams[\"figure.figsize\"] = 10, 5\n",
    "    plt.title(\"Phase portrait of the solution\")\n",
    "    plt.xlabel(\"predator population\")\n",
    "    plt.ylabel(\"prey population\")\n",
    "\n",
    "\n",
    "style = {'description_width': 'initial'}\n",
    "layout = Layout(width=\"400px\")\n",
    "slid1 = widgets.FloatSlider(description='Initial Prey Population $a$ $\\quad$', min=0, max=1, step=0.1, value=0.5,\n",
    "                            style=style, layout=layout)\n",
    "slid2 = widgets.FloatSlider(description='Initial Predator Population $b$', min=0, max=1, step=0.1, value=0.7, style=style,\n",
    "                            layout=layout)\n",
    "interact(plot_pp, a=slid1, b=slid2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [
    {
     "data": {
      "text/plain": "interactive(children=(FloatSlider(value=0.5, description='Initial Prey Population $a$ $\\\\quad$', layout=Layout…",
      "application/vnd.jupyter.widget-view+json": {
       "version_major": 2,
       "version_minor": 0,
       "model_id": "563ae93cbaef42d88fa62f82221d1888"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<function __main__.plot_pp(a, b)>"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def plot_pp(a,b):\n",
    "    m = GEKKO(remote = True)\n",
    "    n = 100\n",
    "    T = 12\n",
    "    m.time = np.linspace(0,T,n)\n",
    "    y0 = 1\n",
    "\n",
    "    x1 = m.Var(value = a)\n",
    "    x2 = m.Var(value = b)\n",
    "    x3 = m.Var(value = (a-y0)**2 + (b-y0)**2)\n",
    "    u = m.Var(value = 0, lb = 0, ub = 1)\n",
    "    p = np.zeros(n)\n",
    "    p[-1] = T\n",
    "\n",
    "    m.Equation(x1.dt()==x1-x1*x2-0.4*x1*u)\n",
    "    m.Equation(x2.dt()==-x2+x1*x2-0.2*x1*u)\n",
    "    m.Equation(x3.dt()==(x1-y0)**2+(x2-y0)**2)\n",
    "\n",
    "    m.Obj(x3)\n",
    "    m.options.IMODE = 6 # optimal control mode\n",
    "\n",
    "    m.solve(disp=False)\n",
    "    plt.subplot(1,2,1)\n",
    "    plt.plot(m.time,x1.value, label = 'Prey')\n",
    "    plt.plot(m.time,x2.value, label = 'Predator')\n",
    "    plt.xlabel(\"time\")\n",
    "    plt.ylabel(\"population\")\n",
    "    plt.title(\"Solution to the control problem with $y_0$ =%.1f\" %y0)\n",
    "    plt.legend()\n",
    "    plt.subplot(1,2,2)\n",
    "    plt.plot(m.time,u.value, color = \"r\")\n",
    "    plt.rcParams[\"figure.figsize\"]=10,5\n",
    "    plt.title(\"Control function values\")\n",
    "    plt.xlabel(\"time\")\n",
    "    plt.ylabel(\"control value\")\n",
    "\n",
    "\n",
    "style = {'description_width': 'initial'}\n",
    "layout = Layout(width = \"400px\")\n",
    "slid1 = widgets.FloatSlider(description = 'Initial Prey Population $a$ $\\qquad$',min=0.4, max=0.8, step=0.05, value=0.5,style = style, layout = layout)\n",
    "slid2 = widgets.FloatSlider(description = 'Initial Predator Population $b$',min=0.4, max=0.8, step=0.05, value=0.5, style = style, layout = layout)\n",
    "interact(plot_pp, a= slid1 ,b=slid2)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}