{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Contrôle de l'orbite d'un satellite\n", "\n", "On assimile le satellite à un point matériel de masse $m$. Sa position $x=(x_1,x_2)$, dans un système de coordonnées plan centré sur la Terre, est donnée par la solution de l'équation de Newton:\n", "\n", "$$\n", "\\ddot x(t) = - \\mu \\frac{x(t)}{|x(t)|^3} + \\frac{F(t)}{m}\n", "$$\n", "\n", "où $F=(F_1,F_2)$ est la poussée et $\\mu$ la constante gravitationnelle de la Terre. \n", "\n", "Nous allons commencer par tracer quelques trajectoires. Nous aurons besoin des modules `numpy` pour le calcul numérique, `matplotlib.pyplot` pour les graphiques, et de la fonction `odeint` du module `scipy.integrate` pour intégrer les systèmes différentiels. Nous définissons également quelques variables globales:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt \n", "from scipy.integrate import odeint\n", "\n", "p0 = 15.625 # Low Earth Orbit altitude\n", "pf = 42.165 # Geostationary Orbit altitude\n", "mu = 398600.47 * 10**(-9) * 3600**2 # Earth gravitation constant\n", "alpha = np.sqrt(mu / pf**3) # rescaled constant\n", "qb = p0/(1.75*pf) # Rescaled Low Earth Orbit altitude\n", "Tfinal = 24\n", "dt = 0.1# time step " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "En adimensionnant le système afin que l'orbite géostationnaire soit à l'altitude 1, l'équation d'état s'écrit:\n", "\n", "$$\n", "\\left\\{\n", "\\begin{array}{rcl}\n", "\\dot q(t) &=& \\eta(t) \\\\\n", "\\dot \\eta(t) & = & -\\alpha^2\\displaystyle{\\frac{q(t)}{|q(t)|^3}}\n", "\\end{array}\n", "\\right .\n", "$$\n", "\n", "où $q=(q_1,q_2)$ est la position et $\\eta=(\\eta_1,\\eta_2)$ la vitesse du satellite. La fonction suivante implémente le champ de vecteur qui définit ce système dynamique :\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def motion_eq(q, t):\n", " \"\"\"\n", " Motion equations for the satellite trajectory.\n", " \n", " Args:\n", " q[0], q[1]: position\n", " q[2], q[3]: velocity\n", "\n", " Returns: \n", " Trajectory equations\n", " \"\"\"\n", " q1,q2,eta1,eta2 = q[0],q[1],q[2],q[3]\n", " k = alpha**2 / np.sqrt( q1**2 + q2**2 )**3\n", " return [eta1, \n", " eta2, \n", " - k * q1, \n", " - k * q2]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Tracé de quelques trajectoires\n", "\n", "On construit une liste contenant les paramètres de quelques trajectoires :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#\n", "test_case = []\n", "# test_case format: (\"label curve\",\"style curve\",q_initial,T_final_plot)\n", "# where q_initial contains the initial position and velocity [q_1, q_2, eta_1, eta_2]\n", "test_case.append((\"Geostationary orbit\",\"b\",[0, 1, -alpha, 0],24))\n", "test_case.append((\"Circular Low Earth Orbit\",\"g\",[qb, 0, 0, alpha/np.sqrt(qb)],3))\n", "test_case.append((\"Elliptic Orbit 1\",\"r\",[-p0/(0.25*pf), 0, 0, -0.2*alpha/np.sqrt(p0/pf)],20))\n", "test_case.append((\"Elliptic Orbit 2\",\"y\",[p0/(1.75*pf), 0, 0, 1.25*alpha/np.sqrt(p0/(1.75*pf))],14))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On parcourt alors la liste des cas tests, on intègre l'équation d'état avec la fonction `odeint` et on trace les trajectoires et les énergies." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, ax = plt.subplots(1, 1, figsize=(10,10))\n", "fig2, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(18,12))\n", "plt.rc('font',size=18)\n", "\n", "\n", "for case in test_case:\n", " case_label,case_style,q0,T_final_plot = case[0],case[1],case[2],case[3]\n", " ts = np.linspace(0, Tfinal, int(Tfinal/dt)+1)\n", "\n", " # Differential equation solution\n", " qs = odeint(motion_eq, q0, ts)\n", " \n", " # Plot of the trajectories\n", " position_x,position_y = qs[:,0],qs[:,1]\n", " vitesse_x,vitesse_y = qs[:,2],qs[:,3]\n", " ax.set_title(\"Various trajectories\")\n", " ax.set_xlabel(\"x\")\n", " ax.set_ylabel(\"y\")\n", " ax.set_aspect('equal')\n", " ax.set_xlim(-1.3,1.5)\n", " ax.set_ylim(-1.2,1.7)\n", " n_plot = (int)(T_final_plot/dt)+1\n", " ax.plot(position_x[0:n_plot], position_y[0:n_plot],case_style+\"+\",label=case_label)\n", " ax.legend()\n", "\n", " # Plot of the energies\n", " Ec = 0.5 * (vitesse_x * vitesse_x + vitesse_y * vitesse_y)\n", " Ep = - alpha**2 / np.sqrt(position_x * position_x + position_y * position_y) \n", " ax1.set_title(\"Energies\")\n", " ax1.set_ylabel(\"Kinetic energy\")\n", " ax1.plot(ts,Ec,case_style+'-')\n", " ax2.set_ylabel(\"Potential energy\")\n", " ax2.plot(ts,Ep,case_style+'-')\n", " ax3.set_ylabel(\"Total energy\")\n", " ax3.plot(ts,Ec+Ep,case_style+'-')\n", " ax3.set_xlabel(\"Time\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Transfert en poussée minimale " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On veut passer d'une orbite basse à une orbite géostationnaire en deux manoeuvres impulsionnelles ([transfert de Hohman](https:\\\\fr.wikipedia.org\\wiki\\Orbite_de_transfert)), c'est-à-dire en imposant deux changements de vitesse, $\\Delta V_0$ à l'instant initial, et $\\Delta V_f$ à l'instant final $t_f$. \n", "\n", "On veut donc que la vitesse initiale soit $\\eta(0) = \\eta_0 + \\Delta V_0$, et qu'au temps $t_f$, on soit sur une orbite géostationnaire, ce qui se traduit par $|q(t_f)|=1$ et $\\eta(t_f) + \\Delta V_f = \\alpha {\\cal R}q(t_f)$, où $\\cal R$ est la rotation d'angle $\\pi/2$. Le but est de minimiser $|\\Delta V_0|^2 + |\\Delta V_f|^2$. \n", "\n", "On cherche donc à résoudre le problème de minimisation:\n", "$$\n", "\\min_{u_0\\in{\\mathbb R}^2} J(u_0) = \\min_{u_0\\in{\\mathbb R}^2} \\frac{1}{2}\\left( |u_0-\\eta_0 |^2 + |\\eta(t_f) - \\alpha{\\cal R}q(t_f)|^2 \\right ),\n", "$$\n", "sous la contrainte que $(q,\\eta$) vérifie les équations d'état sur $[0,t_f]$, avec comme donnée initiale $(q_0,u_0)$, tout en atteignant la **cible** $|q(t_f)|=1$. Noter que l'instant final $t_f$ fait partie des inconnues." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Commençons par tracer les orbites de départ et d'arrivée, qui correspondent aux paramètres suivants: \n", "\n", "|Orbite |Position initiale|  Vitesse initiale         |\n", "|---------------|---------------|----------------------------------------|\n", "|Géostationnaire |$q(0)=(0,1)$ | $\\eta(0)=(-\\alpha,0)$ |\n", "|Circulaire basse |$q(0)=(q_b,0)$ | $\\eta(0)=\\left(0,\\frac{\\alpha}{\\sqrt{q_b}}\\right)$ |\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "q1_0, q2_0 = qb, 0\n", "q1dot_0, q2dot_0= 0, alpha/np.sqrt(qb)\n", "test_case = []\n", "# test_case format: (\"label curve\",\"style curve\",q_initial,T_final_plot)\n", "test_case.append((\"Geostationary orbit\",\"b\",[0, 1, -alpha, 0],24))\n", "test_case.append((\"Circular Low Earth Orbit\",\"g\",[q1_0,q2_0,q1dot_0,q2dot_0],2.5))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(10,10))\n", "plt.rc('font',size=18)\n", "\n", "for case in test_case:\n", " case_label,case_style,q0,T_final_plot = case[0],case[1],case[2],case[3]\n", " ts = np.linspace(0, Tfinal, int(Tfinal/dt)+1)\n", "\n", " # ODE solution\n", " qs = odeint(motion_eq, q0, ts)\n", " \n", " # Plot of the orbit\n", " position_x,position_y = qs[:,0],qs[:,1]\n", " vitesse_x,vitesse_y = qs[:,2],qs[:,3]\n", " ax.set_title(\"A Low Earth Orbit and a Geostationary Orbit\")\n", " ax.set_xlabel(\"x\")\n", " ax.set_ylabel(\"y\")\n", " ax.set_aspect('equal')\n", " ax.set_xlim(-1.3,1.3)\n", " ax.set_ylim(-1.2,1.6)\n", " n_plot = (int)(T_final_plot/dt)+1\n", " ax.plot(position_x[0:n_plot], position_y[0:n_plot],case_style+\"+\",label=case_label)\n", " ax.legend()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### État adjoint\n", "\n", "On introduit l'**état adjoint** $p(t)=(s(t),\\mu(t))\\in{\\mathbb R}^4$ solution de $\\dot p = -A(q,\\eta)^T p$ où\n", "\n", "$$\n", "A = \\begin{bmatrix} \n", "0_{{\\mathbb R}^{2\\times 2}} & I_{{\\mathbb R}^{2\\times 2}} \\\\ \n", "\\Gamma & 0_{{\\mathbb R}^{2\\times 2}}\n", "\\end{bmatrix}\\in {\\mathbb R}^{4\\times 4}\n", "\\qquad \n", "\\Gamma = \\frac{\\alpha^2}{|q|^5}\\begin{bmatrix} \n", "2 q_1^2 - q_2^2 & 3 q_1 q_2 \\\\ \n", " 3 q_1 q_2 & 2 q_2^2 - q_1^2 \n", "\\end{bmatrix}\\in {\\mathbb R}^{2\\times 2}\n", "$$\n", "\n", "La fonction suivante définit le champ de vecteurs de $\\mathbb R^8$ du système d'équations différentielles satisfait par l'état $y(t)=(q(t),\\eta(t))$ et l'état adjoint $p(t) =(s(t),\\mu(t))$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def motion_adjoint_eq_minimum_thrust(q, t):\n", " \"\"\"\n", " Motion equations for the satellite trajectory\n", " and adjoint equations for the minimum thrust problem.\n", " \n", " Args:\n", " q[0], q[1]: position (q)\n", " q[2], q[3]: velocity (eta)\n", " q[4], q[5]: adjoint position (s) \n", " q[6], q[7]: adjoint velocity (mu)\n", " \n", " Returns:\n", " Motion and adjoint equations.\n", " \"\"\"\n", " q1,q2,eta1,eta2 = q[0],q[1],q[2],q[3]\n", " s1,s2,mu1,mu2 = q[4],q[5],q[6],q[7]\n", " tmp = np.sqrt( q1**2 + q2**2 )\n", " k1 = alpha**2 / tmp**3\n", " k2 = alpha**2 / tmp**5\n", " return [eta1, \n", " eta2, \n", " - k1 * q1, \n", " - k1 * q2,\n", " - k2 * ((2*q1**2 - q2**2)*mu1 + (3*q1*q2 )*mu2),\n", " - k2 * ((3*q1*q2 )*mu1 + (2*q2**2 - q1**2)*mu2),\n", " - s1,\n", " - s2\n", " ]\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Méthode de tir\n", "\n", "Le principe d'une **méthode de tir** est de chercher de manière itérative les conditions initiales sur l'état et l'état adjoint qui permettent de vérifier les conditions à l'instant final.\n", "\n", "Ici, la condition initiale sur l'état est donnée par $(q(0),\\eta(0))$, où $q(0)=q_0$ est connue et $\\eta(0) = u_0$ est inconnue (c'est le contrôle). La condition initiale sur l'état adjoint est $(s_0,\\mu_0)$, où $s_0$ est inconnue et $\\mu_0=\\eta_0 - u_0$ d'après la condition d'optimalité. Ainsi la donnée des cinq scalaires $u_0\\in{\\mathbb R}^2$, $s_0\\in{\\mathbb R}^2$ et $t_f$ permet de résoudre le système différentiel vérifié par l'état et l'état adjoint sur $[0,t_f]$. \n", "\n", "À l'instant final $t_f$, les variables d'état $(q,\\eta)$ et adjointes $(s,\\mu)$ vérifient les relations suivantes:\n", "\n", "- Atteinte de la cible :\n", "$$\n", "|q(t_f)|^2 - 1 = 0\n", "$$\n", "\n", "- Conditions de transversalité :\n", "\\begin{align}\n", "-\\mu(t_f) + \\eta(t_f) - \\alpha {\\cal R} q(t_f) & = 0, \\\\\n", "( - s(t_f) + \\alpha {\\cal R}\\eta(t_f) , \\cal R q(t_f))_{{\\mathbb R}^2} & = 0, \\\\\n", "\\end{align}\n", "\n", "- Annulation du hamiltonien :\n", "$$\n", "( s(t_f), \\eta(t_f))_{{\\mathbb R}^2} - \\alpha^2 ( \\mu(t_f), q(t_f))_{{\\mathbb R}^2} = 0.\n", "$$\n", "\n", "La fonction python suivante permet d'évaluer\n", "$$\n", "F_p(u_0,s_0,t_f) \n", "= \n", "\\begin{pmatrix}\n", "|q(t_f)|^2 - 1\\\\\n", "-\\mu(t_f) + \\eta(t_f) - \\alpha {\\cal R} q(t_f) \\\\\n", "( - s(t_f) + \\alpha {\\cal R}\\eta(t_f) , \\cal R q(t_f))_{{\\mathbb R}^2} \\\\\n", "( s(t_f), \\eta(t_f))_{{\\mathbb R}^2} - \\alpha^2 ( \\mu(t_f), q(t_f))_{{\\mathbb R}^2} \\end{pmatrix}\n", "$$\n", "\n", "Noter que chaque évaluation de $F_p$ nécessite une résolution complète du système différentiel état + adjoint. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def Fp_minimum_thrust(X):\n", " \"\"\"\n", " Function to cancel (shooting method) for the minimum thrust problem.\n", " \n", " Args:\n", " X = [eta1_init, eta2_init, s1_init, s2_init, tf]\n", "\n", " Return:\n", " The 5 quantities to cancel.\n", " \"\"\"\n", " global nb_eval,resultats\n", " t1 = np.linspace(0,X[4],int(X[4]/dt)+1)\n", " # The equation is integrated until the final time tf\n", " qs = odeint(motion_adjoint_eq_minimum_thrust, \n", " [q1_0, # initial position (q1)\n", " q2_0, # initial position (q2)\n", " X[0], # initial velocity (eta1 = u0_1 = control)\n", " X[1], # initial velocity (eta2 = u0_2 = control)\n", " X[2], # initial adjoint position (s1)\n", " X[3], # initial adjoint position (s2)\n", " -X[0]+q1dot_0, # initial adjoint velocity (mu1), q1dot_0 = eta0_1 (given original velocity of the trajectory)\n", " -X[1]+q2dot_0],# initial adjoint velocity (mu2), q2dot_0 = eta0_2 (given original velocity of the trajectory)\n", " # given by the optimality condition\n", " t1)\n", " # qs is a np.array with 8 columns and as many rows as time steps \n", " # The last row is extracted to get the solution at time t=tf\n", " q1,q2,eta1,eta2,s1,s2,mu1,mu2 = qs[-1,:] \n", " Rq1,Rq2 = -q2,q1\n", " Reta1,Reta2 = -eta2,eta1\n", " tmp = np.sqrt( q1**2 + q2**2 )\n", " k1 = alpha**2 / tmp**3\n", " XZero = [0,0,0,0,0]\n", " XZero[0] = q1**2 + q2**2 - 1\n", " XZero[1] = -mu1 + eta1 - alpha * Rq1\n", " XZero[2] = -mu2 + eta2 - alpha * Rq2\n", " XZero[3] = (-s1 + alpha*Reta1)*Rq1 + (-s2 + alpha*Reta2)*Rq2\n", " XZero[4] = s1*eta1 + s2*eta2 - k1 * (mu1*q1 + mu2*q2)\n", " \n", " # For the post-processing, the following data are stored:\n", " # iteration number; trajectory (q1,q2); cost value; X\n", " Jval = 0.5 * ( (q1dot_0 - X[0])**2 + (q2dot_0 - X[1])**2 ) + 0.5 * ( (eta1 - alpha*Rq1)**2 + (eta2 - alpha*Rq2)**2 )\n", " nb_eval = nb_eval + 1\n", " resultats.append((nb_eval,qs[:,0],qs[:,1],Jval,X[0],X[1],X[2],X[3],X[4]))\n", " return XZero\n", " \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Le but de la méthode de tir est de chercher $(u_0,s_0,t_f)\\in{\\mathbb R}^5$ tel que\n", "\n", "$$\n", "F_p(u_0,s_0,t_f) = 0\n", "$$\n", "\n", "On va utiliser pour cela la fonction `fsolve`, du module `scipy.optimize`, qui permet de trouver un zéro d'une fonction." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import fsolve" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On pourra vérifier que l'initialisation de la méthode de tir de tir est très délicate. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "u0_1, u0_2 = -0.6, 0.6\n", "s0_1, s0_2 = 1., 1.\n", "tf = 1.\n", "nb_eval = 0 # how many times the state and adjoint equations are solved (postprocessing)\n", "resultats = [] # to store the intermediate results at each shooting iteration (postprocessing)\n", "\n", "a = fsolve(Fp_minimum_thrust,[u0_1,u0_2,s0_1,s0_2,tf])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Affichage des résultats\n", "\n", "On trace l'évolution de la fonction coût $J$, et des 5 inconnues $\\eta_1(0) = u_{0,1},\\eta_2(0)= u_{0,2},s_1(0),s_2(0),t_f$ au cours des itérations de la méthode de tir:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot of the cost function and shooting variables\n", "fig1,(ax1,ax2,ax3) = plt.subplots(3, 1,figsize=(16, 10))\n", "fig2,(ax4,ax5,ax6) = plt.subplots(3, 1,figsize=(16, 10))\n", "plt.rc('font',size=18)\n", "\n", "Jval = []\n", "eta1val = []\n", "eta2val = []\n", "s1val = []\n", "s2val = []\n", "tfval = []\n", "\n", "for iter,position_x,position_y,Jv,eta1,eta2,s1,s2,tf in resultats:\n", " Jval.append(Jv)\n", " eta1val.append(eta1)\n", " eta2val.append(eta2)\n", " s1val.append(s1)\n", " s2val.append(s2)\n", " tfval.append(tf)\n", "\n", "ax1.set_ylabel(\"Cost function\")\n", "ax1.plot(Jval,'r')\n", "\n", "ax2.set_ylabel(\"Final time\")\n", "ax2.plot(tfval,'g')\n", "\n", "ax3.set_ylabel(\"eta1\")\n", "ax3.plot(eta1val,'b')\n", "\n", "\n", "ax4.set_ylabel(\"eta2\")\n", "ax4.plot(eta2val,'b')\n", "\n", "ax5.set_ylabel(\"s1\")\n", "ax5.plot(s1val,'y')\n", "\n", "ax6.set_ylabel(\"s2\")\n", "ax6.plot(s2val,'y')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On trace la trajectoire initiale et celle obtenue à convergence:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=[12,12])\n", "plt.rc('font',size=18)\n", "#\n", "# Plot of the two orbits (initial and target)\n", "#\n", "for case in test_case:\n", " case_label,case_style,q0,T_final_plot = case[0],case[1],case[2],case[3]\n", " ts = np.linspace(0, Tfinal, int(Tfinal/dt)+1)\n", "\n", " # ODE solution \n", " qs = odeint(motion_eq, q0, ts)\n", " \n", " # Plot of the orbit\n", " position_x,position_y = qs[:,0],qs[:,1]\n", " vitesse_x,vitesse_y = qs[:,2],qs[:,3]\n", " ax.set_xlabel(\"x\")\n", " ax.set_ylabel(\"y\")\n", " ax.set_aspect('equal')\n", " ax.set_xlim(-1.3,1.3)\n", " ax.set_ylim(-1.2,1.5)\n", " n_plot = (int)(T_final_plot/dt)+1\n", " ax.plot(position_x[0:n_plot], position_y[0:n_plot],case_style+\"-\",label=case_label)\n", "\n", "# Plot of the trajectory corresponding to the initial parameters given to the shooting algorithm\n", "\n", "iter,position_x,position_y,Jv,eta1,eta2,s1,s2,tf = resultats[0]\n", "ax.plot(position_x,position_y,'orange',label='Initial guess')\n", "\n", "# Plot of the optimal trajectory (last iteration of the shooting algorithm)\n", "\n", "iter,position_x,position_y,Jv,eta1,eta2,s1,s2,tf = resultats[-1]\n", "ax.plot(position_x,position_y,'red',label='Optimal trajectory')\n", "\n", "ax.set_title(\"Optimal thrust trajectory from a Low Earth to a Geostationary Orbit\")\n", "ax.legend()\n", "\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Animation des trajectoires successives obtenues par l'algorithme de tir. **N.B.:** nécessite que la librairie `ffmpeg` soit installée sur votre ordinateur." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib import animation, rc\n", "from IPython.display import HTML\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=[12,12])\n", "plt.rc('font',size=18)\n", "#\n", "# Plot of the two orbits (initial and target)\n", "#\n", "for case in test_case:\n", " case_label,case_style,q0,T_final_plot = case[0],case[1],case[2],case[3]\n", " ts = np.linspace(0, Tfinal, int(Tfinal/dt)+1)\n", "\n", " # ODE solution \n", " qs = odeint(motion_eq, q0, ts)\n", " \n", " # Plot of the orbit\n", " position_x,position_y = qs[:,0],qs[:,1]\n", " vitesse_x,vitesse_y = qs[:,2],qs[:,3]\n", " ax.set_xlabel(\"x\")\n", " ax.set_ylabel(\"y\")\n", " ax.set_aspect('equal')\n", " ax.set_xlim(-1.3,1.3)\n", " ax.set_ylim(-1.2,1.5)\n", " n_plot = (int)(T_final_plot/dt)+1\n", " ax.plot(position_x[0:n_plot], position_y[0:n_plot],case_style+\"-\",label=case_label)\n", "\n", "\n", "line, = ax.plot([], [], lw=2)\n", "\n", "def init():\n", " line.set_data([], [])\n", " return (line,)\n", "\n", "def animate(i):\n", " x = resultats[i-1][1]\n", " y = resultats[i-1][2]\n", " line.set_data(x, y)\n", " return (line,)\n", "\n", "# call the animator. blit=True means only re-draw the parts that have changed.\n", "anim = animation.FuncAnimation(fig, animate, init_func=init,\n", " frames=100, interval=20, blit=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "HTML(anim.to_jshtml()) # if this does not work try HTML(anim.to_html5_video())" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }