{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "54086665-879b-438f-af50-2ee8ccb39d49",
   "metadata": {},
   "source": [
    "## MT09 - TP5 - Automne 2024\n",
    "### Approximation aux moindres carrés\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5dd0ce4-cb08-4aba-a9ac-bd51b794385e",
   "metadata": {},
   "source": [
    "#### Ex 1 - Approximation affine par morceaux\n",
    "\n",
    "On considère le nuage de $m$ points $(x_i,y_i)_{i=1,...,m}$ suivants donnés par les tableaux $x$ et $y$ du script code python ci-dessous :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d688136-b21c-44fb-a51f-a99f4f62a3f6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "m = 20\n",
    "x = np.linspace(0.0, 5.0, m) + 0.05*np.random.randn(m)\n",
    "y = 0.8 + 0.5*x-x**2/6 + 0.2*np.random.randn(m)\n",
    "plt.plot(x, y, 'o'); plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54b3f1d5-16ba-4ade-8fb9-787fc1146231",
   "metadata": {},
   "source": [
    "On souhaite représenter ce nuage de points par une ligne brisée en minimisant les écarts aux carrés entre les points et la ligne brisée.\n",
    "Dans un premier temps, coder la fonction dite chapeau ```y = phi(x)``` définie par\n",
    "\n",
    "$$\n",
    "\\phi(x) = \\max(0, 1-|x|).\n",
    "$$\n",
    "\n",
    "Représenter graphiquement la fonction $\\phi$ sur l'intervalle $[-3,3]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d356719-f874-4991-ab99-db793ff66ead",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "95453e70-1bf1-48a9-8b93-8d9ce34e40d8",
   "metadata": {},
   "source": [
    "On cherche ensuite une approximation $f(x)$ de la forme\n",
    "\n",
    "$$\n",
    "f(x) = \\sum_{j=0}^5 \\, u_j\\, \\phi(x - j),\n",
    "$$\n",
    "où les $(u_j)_{j=0,...,5}$ sont les inconnues. Réécrire le problème aux moindres carrés\n",
    "\n",
    "$$\n",
    "\\min_{(u_0,...,u_5)} \\quad \\frac{1}{2}\\sum_{i=1}^m\\, [f(x_i)-y_i]^2\n",
    "$$\n",
    "\n",
    "sous la forme\n",
    "\n",
    "$$\n",
    "\\min_{u\\in\\mathbb{R}^6} \\quad \\frac{1}{2} \\|Au - y\\|^2\n",
    "$$\n",
    "en déterminant la matrice $A$. Coder en python la matrice $A$, résoudre les équations normales et tracez sur le même graphique le nuage de points et la fonction d'approximation $f(x)$ ainsi trouvée. Afficher le vecteur $\\mathbf{u}$. Comment interprétez-vous les valeurs $u_i$ de $\\mathbf{u}$ ?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a4618a0-88c6-4cd5-be76-909561690c97",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "833c28d6-c8ec-4589-af4a-12c87100048e",
   "metadata": {},
   "source": [
    "#### Ex2 - Recherche d'une ellipse passant par un nuage de points\n",
    "\n",
    "On considère le nuage de points $(x_i,y_i)$, $i=1,...,m$ générés par le script python ci-dessous :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d84359f-80e0-4bcc-9168-7ff9880dee37",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "a = 2.5\n",
    "b = 1.1\n",
    "m = 50\n",
    "theta = np.linspace(0, 2*np.pi, m)\n",
    "xi = a * np.cos(theta) + 0.05*np.random.randn(m)\n",
    "yi = b * np.sin(theta) + 0.05*np.random.randn(m)\n",
    "plt.plot(xi, yi, 'or'); plt.axis(\"equal\"); plt.grid()\n",
    "plt.xlim([-3, 3]); plt.xlabel('x'); plt.ylabel('y');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb8405e5-978b-4c48-b76e-06857d5e8b32",
   "metadata": {},
   "source": [
    "Par ces points, on va essayer de faire passer \"au mieux' une ellipse de paramètres d'axe\n",
    "$\\hat a$ et $\\hat b$, d'équation \n",
    "\n",
    "$$\n",
    "\\frac{x^2}{\\hat a^2} + \\frac{y^2}{\\hat b^2} = 1\n",
    "$$\n",
    "\n",
    "ou, de manière équivalente, d'équations paramétrées\n",
    "\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "x(\\theta) = \\hat a \\cos(\\theta), \\\\\n",
    "y(\\theta) = \\hat b \\sin(\\theta), \\quad \\theta\\in[0, 2\\pi].\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "On pourrait chercher $(\\hat a,\\hat b)$ qui réalisent le minimum de\n",
    "\n",
    "$$\n",
    "\\min_{(a,b)} \\quad \\frac{1}{2} \\sum_{i=1}^m \\left[ \\frac{x_i^2}{\\hat a^2} + \\frac{y_i^2}{\\hat b^2}-1\\right]^2\n",
    "$$\n",
    "\n",
    "mais ce problème ne peut pas être mis sous la forme $\\min_u \\frac{1}{2} \\|Au-b\\|^2$. On considère plutôt le problème aux moindres carrés\n",
    "\n",
    "$$\n",
    "\\min_{u=(u_1,u_2)} \\quad \\frac{1}{2} \\sum_{i=1}^m \\left[ x_i^2\\, u_1 + y_i^2\\, u_2 -1\\right]^2\n",
    "$$\n",
    "puis on affectera $\\hat a = \\frac{1}{\\sqrt{u_1}}$ et $\\hat b = \\frac{1}{\\sqrt{u_2}}$. Ecrire ce problème de minimisation sous la forme\n",
    "\n",
    "$$\n",
    "\\min_u \\quad \\frac{1}{2}\\, \\|Au-b\\|^2\n",
    "$$\n",
    "\n",
    "où l'on déterminera la matrice $A$ et le second membre $b$. Résoudre les équations normales, calculer\n",
    "$\\hat a$ et $\\hat b$. Sur le même graphique, tracer le nuage de points $(x_i,y_i)$ et l'ellipse de pramètres\n",
    "d'axe $\\hat a$ et $\\hat b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da6e3f75-6ced-432b-84af-f7dd03d90b9a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "4eebd9da-9b86-40fd-8ab7-26ec04392ca4",
   "metadata": {},
   "source": [
    "#### Ex 3 - Identification de paramètres d'un système masse-ressort\n",
    "\n",
    "On considère ici un système masse-ressort linéaire avec amortissement, $m$ est la masse, $k$ est le coefficient de raideur du ressort et $a$ est le coefficient d'amortissement. Le mouvement du point matériel est gouverné par l'équation différentielle du second ordre:\n",
    "\n",
    "$$\n",
    "m\\, \\ddot{x}(t) + a\\, \\dot x(t) + k\\, x(t) = 0\n",
    "$$\n",
    "\n",
    "où $x(t)$ représente la position de la masse au temps $t$.\n",
    "L'équation caractéristique correspondante est:\n",
    "\n",
    "$$\n",
    "r^2 + \\frac{a}{m} r + \\frac{k}{m} = 0\n",
    "$$\n",
    "\n",
    "avec le discriminant $\\Delta$ fonction de $m$, $k$ et $a$ :\n",
    "\n",
    "$$\n",
    "\\Delta = \\left(\\frac{a}{m}\\right)^2 - 4 \\frac{k}{m}.\n",
    "$$\n",
    "\n",
    "On suppose $\\Delta<0$, si bien que la masse à un mouvement oscillatoire amorti. Sans conditions initiales, la solution générale est de la forme\n",
    "\n",
    "$$\n",
    "x(t) = A\\, \\exp(-\\frac{a}{2m}t)\\, \\cos(\\frac{\\sqrt{-\\Delta}}{2} t) \n",
    "+ B\\, \\exp(-\\frac{a}{2m}t)\\, \\sin(\\frac{\\sqrt{-\\Delta}}{2} t) \n",
    "$$\n",
    "\n",
    "avec deux constantes réelles arbitraires $A$ et $B$. \n",
    "\n",
    "Dans la suite, on considère $m=1$, $a=0.4$ et $k=11$ et on suppose $A=1$ et $B=0$ correspondant à $x(0)=1$ et $\\dot x(0) = -0.2$. La solution est donc\n",
    "\n",
    "$$\n",
    "x(t) = \\exp(-\\frac{a}{2m}t)\\, \\cos(\\frac{\\sqrt{-\\Delta}}{2} t).\n",
    "$$\n",
    "\n",
    "Sur brouillon, calculez $\\dot x(t)$ et $\\ddot x(t)$. Ecrire une fonction \n",
    "\n",
    "```\n",
    "x, xdot, xdotdot = xressort(t)\n",
    "``` \n",
    "\n",
    "qui retourne les valeurs $x(t)$, $\\dot x(t)$ et $\\ddot x(t)$ au temps $t$. Afficher le résultat de\n",
    "```xressort(0)```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23051076-08f4-4714-a60c-6af738f1633a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b573f014-e3cd-4dfb-97fd-284cc0dd5a3f",
   "metadata": {},
   "source": [
    "L'exercice consiste à étudier si, à partir de la connaissance de la position $x(t^n)$, la vitesse $\\dot x(t^n)$ et l'accélération $\\ddot x(t^n)$ à différents instants $t^n$, **connaissant la masse $m$**, on est en mesure d'identifier les caractéristiques mécaniques du ressort, c'est-à-dire la **raideur $k$** et l'**amortissement $a$**. Pour chaque instant $t^n$, on a donc\n",
    "\n",
    "$$\n",
    "\\ddot x(t^n) = - \\frac{a}{m}\\, \\dot x(t^n) - \\frac{k}{m}\\, x(t^n)\n",
    "$$\n",
    "\n",
    "Les mesures de $x(t)$, $\\dot x(t)$ et $\\ddot x(t)$ ont sujettes à des incertitudes de mesures représentées par un bruit gaussien. Considérez le code python ci-dessous pour la génération de données. Dans ce cas,\n",
    "\n",
    "$$\n",
    "\\ddot x(t^n) \\approx - \\frac{a}{m}\\, \\dot x(t^n) - \\frac{k}{m}\\, x(t^n).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1749152-3a5a-443b-8f95-785f6f147095",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "N = 200\n",
    "t = np.linspace(0, 10, N)\n",
    "x, xdot, xdotdot = xressort(t)\n",
    "x = x*(1 + 0.02*np.random.randn(N))\n",
    "xdot = xdot*(1 + 0.02*np.random.randn(N))\n",
    "xdotdot = xdotdot*(1 + 0.02*np.random.randn(N))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe994d8f-e50a-42d1-b684-30f33b7d5374",
   "metadata": {},
   "source": [
    "On considère le problème aux moindres carrés\n",
    "\n",
    "$$\n",
    "\\min_{u=(k,a)}\\quad \\frac{1}{2} \\sum_{n=1}^N \\left[ \\ddot x(t^n) + \\frac{a}{m} \\dot x(t^n)\n",
    "+ \\frac{k}{m} x(t^n)\\right]^2.\n",
    "$$\n",
    "\n",
    "Ecrire le problème sous la forme matricielle\n",
    "\n",
    "$$\n",
    "\\min_{u=(k,a)} \\quad \\frac{1}{2} \\|A u - b\\|^2\n",
    "$$\n",
    "\n",
    "en déterminant la matrice $A$ et le second membre $b$. Résoudre les équations normales et affichez les valeurs\n",
    "estimées $(\\hat k, \\hat a)$ et $k$ et $a$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d74949b-3850-4d1a-9848-4e66ae4aaf22",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "2eb21df5-e152-4777-9ffc-d251cc544ef7",
   "metadata": {},
   "source": [
    "Sur le même graphique, tracer la solution exacte $x(t)$ et $\\hat x(t)$ correspondant à la solution de paramètres mécaniques estimés $\\hat k$ et $\\hat a$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d0ddf4b-9ca2-4e96-9ecb-7c92bd24db97",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "528e7d74-76d4-4192-adc8-9805891a7647",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}