{ "cells": [ { "cell_type": "markdown", "id": "190011e1", "metadata": {}, "source": [ "# Simulación de Monte Carlo de un sistema estocástico: confiabilidad de un sistema de 5 componentes\n", "\n", "**Objetivo.** Modelar el tiempo de falla de un sistema de 5 componentes (Law & Kelton, Ejemplo 13.9) a partir de su diagrama de confiabilidad, estimar $E[Y]$ y $P(Y \\ge 7)$ mediante simulación de Monte Carlo, y visualizar dinámicamente tanto el diagrama de bloques como la evolución del sistema en el tiempo.\n", "\n", "Curso: Modelado de Sistemas bajo Incertidumbre · Universidad de los Andes · 2026-20\n" ] }, { "cell_type": "code", "execution_count": null, "id": "da2fc0fc", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from matplotlib.patches import Rectangle\n", "from matplotlib import animation\n", "from IPython.display import HTML\n", "\n", "np.random.seed(2026)\n", "\n", "betas = np.array([10.0, 8.0, 7.0, 5.0, 6.0]) # beta_1, ..., beta_5 (medias, en días)\n", "n = 25_000 # número de réplicas\n", "T = 7.0 # umbral para P(Y >= T)\n" ] }, { "cell_type": "markdown", "id": "50fe7b68", "metadata": {}, "source": [ "## 0. Enunciado del ejercicio (Law & Kelton, Ejemplo 13.9)\n", "\n", "Se desea enviar una señal de $A$ a $B$ a través de un sistema con 5 componentes, cada uno sujeto a fallas aleatorias e independientes. Sea $Y$ el tiempo de falla de todo el sistema, y sea $X_i$ el tiempo de falla del componente $i$, para $i=1,\\dots,5$.\n", "\n", "El diagrama de confiabilidad (Fig. 13.14) es: el componente **1** en paralelo con la rama formada por el componente **2** en serie con **3 y 4 en paralelo**; todo eso en serie con el componente **5**.\n", "\n", "| $i$ | $\\beta_i$ (días) |\n", "|---|---:|\n", "| 1 | 10 |\n", "| 2 | 8 |\n", "| 3 | 7 |\n", "| 4 | 5 |\n", "| 5 | 6 |\n", "\n", "Cada $X_i$ se distribuye exponencial con media $\\beta_i$. Puede demostrarse (ver Prob. 13.5 del libro) que\n", "\n", "$$Y = \\min\\Big(\\max\\{X_1,\\ \\min[X_2,\\ \\max(X_3,X_4)]\\},\\ X_5\\Big) \\tag{13.5}$$\n", "\n", "Queremos estimar por simulación: (a) $E[Y]$, el tiempo esperado de falla del sistema, y (b) $p_7 = P(Y \\ge 7)$, la probabilidad de que el sistema siga funcionando después de 7 días. El libro reporta, con $n=25\\,000$ réplicas, $\\bar Y(25\\,000) = 4.33$ días y $\\hat p_7 = 0.19$." ] }, { "cell_type": "markdown", "id": "e2d48f1e", "metadata": {}, "source": [ "## 1. Diagrama del sistema (dado)\n", "\n", "Esta función dibuja el diagrama de confiabilidad de la Fig. 13.14. La vas a reutilizar más adelante para la animación — no hace falta modificarla." ] }, { "cell_type": "code", "execution_count": null, "id": "44b896f3", "metadata": {}, "outputs": [], "source": [ "def caja(ax, x, y, w, h, texto, color):\n", " ax.add_patch(Rectangle((x, y), w, h, facecolor=color, edgecolor=\"black\", lw=1.2, zorder=3))\n", " ax.text(x + w / 2, y + h / 2, texto, ha=\"center\", va=\"center\", fontsize=11, fontweight=\"bold\", zorder=4)\n", "\n", "def linea(ax, x1, y1, x2, y2, color=\"black\", lw=1.4):\n", " ax.plot([x1, x2], [y1, y2], color=color, lw=lw, zorder=2)\n", "\n", "def dibujar_diagrama(ax, colores=None):\n", " \"\"\"Dibuja el diagrama de confiabilidad de la Fig. 13.14. `colores` es un\n", " diccionario opcional {1: color, 2: color, ..., 5: color} para colorear\n", " cada componente (por defecto, todos en azul claro).\"\"\"\n", " if colores is None:\n", " colores = {i: \"#AEC6E8\" for i in range(1, 6)}\n", " ax.clear()\n", " ax.set_xlim(-0.5, 9.5)\n", " ax.set_ylim(-0.5, 5.5)\n", " ax.axis(\"off\")\n", " ax.set_aspect(\"equal\")\n", "\n", " ax.scatter([0], [3], s=250, facecolor=\"white\", edgecolor=\"black\", zorder=5)\n", " ax.text(0, 3, \"A\", ha=\"center\", va=\"center\", fontsize=11, fontweight=\"bold\", zorder=6)\n", " ax.scatter([9], [3], s=250, facecolor=\"white\", edgecolor=\"black\", zorder=5)\n", " ax.text(9, 3, \"B\", ha=\"center\", va=\"center\", fontsize=11, fontweight=\"bold\", zorder=6)\n", "\n", " linea(ax, 0.25, 3, 1, 3)\n", " linea(ax, 1, 3, 1, 4.5)\n", " linea(ax, 1, 3, 1, 1.5)\n", "\n", " linea(ax, 1, 4.5, 2, 4.5)\n", " caja(ax, 2, 4.1, 1.4, 0.8, \"1\", colores[1])\n", " linea(ax, 3.4, 4.5, 7, 4.5)\n", " linea(ax, 7, 4.5, 7, 3)\n", "\n", " linea(ax, 1, 1.5, 2, 1.5)\n", " caja(ax, 2, 1.1, 1.4, 0.8, \"2\", colores[2])\n", " linea(ax, 3.4, 1.5, 4, 1.5)\n", " linea(ax, 4, 1.5, 4, 2.3)\n", " linea(ax, 4, 1.5, 4, 0.7)\n", "\n", " caja(ax, 4, 1.9, 1.4, 0.8, \"3\", colores[3])\n", " caja(ax, 4, 0.3, 1.4, 0.8, \"4\", colores[4])\n", "\n", " linea(ax, 5.4, 2.3, 6, 2.3)\n", " linea(ax, 5.4, 0.7, 6, 0.7)\n", " linea(ax, 6, 2.3, 6, 0.7)\n", " linea(ax, 6, 1.5, 7, 1.5)\n", " linea(ax, 7, 1.5, 7, 3)\n", "\n", " linea(ax, 7, 3, 7.3, 3)\n", " caja(ax, 7.3, 2.6, 1.4, 0.8, \"5\", colores[5])\n", " linea(ax, 8.7, 3, 8.75, 3)\n", "\n", "fig, ax = plt.subplots(figsize=(11, 6))\n", "dibujar_diagrama(ax)\n", "fig.suptitle(\"Diagrama de confiabilidad del sistema (Fig. 13.14)\")\n", "fig.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "d2dc70b0", "metadata": {}, "source": [ "## 2. Pseudocódigo\n", "\n", "```\n", "Entrada: β = (β1, ..., β5) medias de las exponenciales, n réplicas, umbral T\n", "Salida: Ȳ(n) ≈ E[Y], p̂_T ≈ P(Y ≥ T)\n", "\n", "Para j = 1, ..., n:\n", " Para i = 1, ..., 5:\n", " generar Ui ~ U(0,1)\n", " Xi ← -βi · ln(Ui) # tiempo de falla del componente i\n", "\n", " # combinar los tiempos de falla según el diagrama de confiabilidad:\n", " rama_34 ← max(X3, X4) # 3 y 4 en paralelo (redundantes)\n", " rama_2_34 ← min(X2, rama_34) # 2 en serie con (3∥4)\n", " rama_paralela ← max(X1, rama_2_34) # rama 1 en paralelo con [2-(3∥4)]\n", " Yj ← min(rama_paralela, X5) # todo en serie con el componente 5\n", "\n", "Fin Para\n", "\n", "Ȳ(n) ← (1/n) Σ Yj\n", "p̂_T ← (1/n) Σ I(Yj ≥ T)\n", "\n", "Graficar histograma de {Y1, ..., Yn}\n", "```" ] }, { "cell_type": "markdown", "id": "c03f1d02", "metadata": {}, "source": [ "## 3. Generación de tiempos de falla \n", "\n", "Implementa la generación vectorizada de los tiempos de falla $X_i = -\\beta_i \\ln(U_i)$ para los 5 componentes y las `n` réplicas al mismo tiempo." ] }, { "cell_type": "code", "execution_count": null, "id": "5e2a2616", "metadata": {}, "outputs": [], "source": [ "def generar_tiempos_falla(betas, n):\n", " \"\"\"\n", " Genera una matriz de tiempos de falla de forma (n, 5): cada fila es una\n", " réplica j, cada columna i es X_i ~ Exponencial(beta_i), usando\n", " X_i = -beta_i * ln(U_i) con U_i ~ U(0,1).\n", " \"\"\"\n", " # TODO 1: genera una matriz U de forma (n, len(betas)) con np.random.random\n", " # TODO 2: calcula X = -betas * np.log(U) (usa broadcasting: betas tiene forma (5,))\n", " # TODO 3: retorna X\n", " raise NotImplementedError(\"Implementa generar_tiempos_falla\")\n", "\n", "# --- prueba rápida (no la modifiques) ---\n", "X_test = generar_tiempos_falla(betas, 1000)\n", "assert X_test.shape == (1000, 5)\n", "assert np.all(X_test > 0)\n", "print(\"Medias muestrales (deberían acercarse a beta):\", np.round(X_test.mean(axis=0), 2))\n", "print(\"beta: \", betas)\n", "print(\"OK: generar_tiempos_falla pasa la prueba básica.\")\n" ] }, { "cell_type": "markdown", "id": "7409db21", "metadata": {}, "source": [ "## 4. Cálculo del tiempo de falla del sistema\n", "\n", "Implementa la ecuación (13.5) de forma vectorizada, siguiendo el pseudocódigo de la sección 2. La función debe funcionar tanto si `X1,...,X5` son arreglos de longitud `n` (simulación completa) como si son 5 números sueltos (una sola réplica, útil más adelante para la animación)." ] }, { "cell_type": "code", "execution_count": null, "id": "1db636b8", "metadata": {}, "outputs": [], "source": [ "def calcular_Y(X1, X2, X3, X4, X5):\n", " \"\"\"Aplica la ecuación (13.5): Y = min( max(X1, min(X2, max(X3,X4))), X5 ).\"\"\"\n", " # TODO 1: rama_34 = max(X3, X4) -> usa np.maximum\n", " # TODO 2: rama_2_34 = min(X2, rama_34) -> usa np.minimum\n", " # TODO 3: rama_paralela = max(X1, rama_2_34) -> usa np.maximum\n", " # TODO 4: Y = min(rama_paralela, X5) -> usa np.minimum\n", " # TODO 5: retorna Y\n", " raise NotImplementedError(\"Implementa calcular_Y\")\n", "\n", "# --- prueba rápida (no la modifiques) ---\n", "# Caso construido a mano: rama_34=max(3,7)=7; rama_2_34=min(5,7)=5; rama_paralela=max(2,5)=5; Y=min(5,9)=5\n", "Y_test = calcular_Y(2, 5, 3, 7, 9)\n", "assert np.isclose(Y_test, 5)\n", "print(\"OK: calcular_Y pasa la prueba básica. Y_test =\", Y_test)\n" ] }, { "cell_type": "markdown", "id": "365b7dc4", "metadata": {}, "source": [ "## 5. Simulación completa: estimación de E[Y] \n", "\n", "Usa las dos funciones anteriores para generar `n = 25 000` réplicas y estimar $E[Y]$ con $\\bar Y(n)$. Compara contra el valor del libro (4.33 días)." ] }, { "cell_type": "code", "execution_count": null, "id": "f08d5421", "metadata": {}, "outputs": [], "source": [ "# TODO: genera X = generar_tiempos_falla(betas, n)\n", "# TODO: calcula Y = calcular_Y(X[:, 0], X[:, 1], X[:, 2], X[:, 3], X[:, 4])\n", "\n", "\n", "# TODO: calcula Ybar = Y.mean() e imprímelo junto con el valor del libro (4.33 días)\n" ] }, { "cell_type": "markdown", "id": "d90fc144", "metadata": {}, "source": [ "## 6. Estimación de $p_7 = P(Y \\ge 7)$ \n", "\n", "Usa la función indicadora $I_j(7,\\infty) = 1$ si $Y_j \\ge 7$, y estima $p_7$ como el promedio de esos indicadores. Compara contra el valor del libro (0.19)." ] }, { "cell_type": "code", "execution_count": null, "id": "65add57f", "metadata": {}, "outputs": [], "source": [ "# TODO: crea el arreglo indicador = (Y >= T), convertido a float\n", "# TODO: calcula p7_hat = indicador.mean()\n", "# TODO: imprime p7_hat junto con el valor del libro (0.19)\n" ] }, { "cell_type": "markdown", "id": "9f324fd7", "metadata": {}, "source": [ "## 7. Histograma de Y \n", "\n", "Reproduce la Fig. 13.15 del libro: un histograma de las 25 000 $Y_j$ generadas." ] }, { "cell_type": "code", "execution_count": null, "id": "2f75d50a", "metadata": {}, "outputs": [], "source": [ "# TODO: arma un histograma de Y (bins=50)\n", "# TODO: agrega una línea vertical en Ybar y otra en T (umbral) para referencia\n", "# TODO: etiqueta los ejes y el título, agrega una leyenda\n" ] }, { "cell_type": "markdown", "id": "da842543", "metadata": {}, "source": [ "## 8. El sistema a través del tiempo \n", "\n", "Hasta ahora solo calculamos el número $Y$ (el instante en que falla el sistema). Para poder **animar** el diagrama, necesitamos una versión de la misma lógica que responda, para cualquier instante $t$: *¿sigue conectado el sistema de $A$ a $B$ en el instante $t$?*\n", "\n", "Un componente $i$ está funcionando en el instante $t$ si $t < X_i$. Traduce la ecuación (13.5) a lógica booleana, con la misma estructura de ramas del pseudocódigo:\n", "\n", "```\n", "arriba_1 ← (t < X1)\n", "arriba_34 ← (t < X3) OR (t < X4) # 3 y 4 en paralelo\n", "arriba_2_34 ← (t < X2) AND arriba_34 # 2 en serie con (3∥4)\n", "arriba_paralela ← arriba_1 OR arriba_2_34 # rama 1 en paralelo con [2-(3∥4)]\n", "arriba_5 ← (t < X5)\n", "sistema_arriba ← arriba_paralela AND arriba_5\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "3f0e6d8e", "metadata": {}, "outputs": [], "source": [ "def sistema_funciona(t, X1, X2, X3, X4, X5):\n", " \"\"\"Retorna True si el sistema sigue conectado de A a B en el instante t.\"\"\"\n", " # TODO 1: arriba_1 = t < X1\n", " # TODO 2: arriba_34 = (t < X3) or (t < X4)\n", " # TODO 3: arriba_2_34 = (t < X2) and arriba_34\n", " # TODO 4: arriba_paralela = arriba_1 or arriba_2_34\n", " # TODO 5: arriba_5 = t < X5\n", " # TODO 6: retorna arriba_paralela and arriba_5\n", " raise NotImplementedError(\"Implementa sistema_funciona\")\n", "\n", "# --- prueba rápida (no la modifiques) ---\n", "X_prueba = (2.0, 5.0, 3.0, 7.0, 9.0) # el mismo caso de la sección 4: Y = 5\n", "Y_prueba = calcular_Y(*X_prueba)\n", "assert sistema_funciona(Y_prueba - 0.001, *X_prueba) == True\n", "assert sistema_funciona(Y_prueba + 0.001, *X_prueba) == False\n", "print(\"OK: sistema_funciona es consistente con calcular_Y. Y_prueba =\", Y_prueba)\n" ] }, { "cell_type": "markdown", "id": "ab0f955c", "metadata": {}, "source": [ "## 9. Visualización dinámica — Animación del diagrama en el tiempo (dado)\n", "\n", "Con `generar_tiempos_falla`, `calcular_Y` y `sistema_funciona` ya implementadas, podemos animar una réplica concreta: los componentes se van poniendo en rojo a medida que fallan, y el título indica si el sistema sigue **conectado** o ya está **desconectado**. Elegimos a propósito una réplica (semilla 1) que ilustra bien la redundancia: el componente 2 falla primero, pero el sistema sigue funcionando gracias al componente 1, hasta que este también falla." ] }, { "cell_type": "code", "execution_count": null, "id": "1a4c3382", "metadata": {}, "outputs": [], "source": [ "np.random.seed(1)\n", "X_ejemplo = generar_tiempos_falla(betas, 1)[0]\n", "np.random.seed(2026) # restaura la semilla del curso\n", "\n", "Y_ejemplo = calcular_Y(*X_ejemplo)\n", "print(\"Tiempos de falla de esta réplica:\", np.round(X_ejemplo, 2))\n", "print(f\"Y de esta réplica: {Y_ejemplo:.2f} días\")\n", "\n", "VERDE, ROJO = \"#6FCF97\", \"#EB5757\"\n", "t_max = Y_ejemplo * 1.3\n", "frames_t = np.linspace(0, t_max, 50)\n", "\n", "fig, ax = plt.subplots(figsize=(11, 6))\n", "titulo = fig.suptitle(\"\")\n", "\n", "def actualizar(i):\n", " t = frames_t[i]\n", " colores = {k + 1: (VERDE if t < X_ejemplo[k] else ROJO) for k in range(5)}\n", " dibujar_diagrama(ax, colores)\n", " arriba = sistema_funciona(t, *X_ejemplo)\n", " estado = \"CONECTADO (A → B)\" if arriba else \"DESCONECTADO — sistema fallado\"\n", " color_estado = \"#1a7a3c\" if arriba else \"#a4222c\"\n", " titulo.set_text(f\"t = {t:5.2f} días Estado del sistema: {estado}\")\n", " titulo.set_color(color_estado)\n", " return []\n", "\n", "ani = animation.FuncAnimation(fig, actualizar, frames=len(frames_t), blit=False)\n", "plt.close(fig)\n", "HTML(ani.to_jshtml(fps=6))\n" ] }, { "cell_type": "markdown", "id": "f48b2436", "metadata": {}, "source": [ "## 11. Preguntas de discusión\n", "\n", "1. En la réplica animada de la sección 9, el componente 4 falla antes que el componente 1, pero eso no cambia el momento en que el sistema completo falla. ¿Por qué?\n", "2. ¿Qué papel juega específicamente la redundancia de los componentes 3 y 4 en la confiabilidad del sistema? Si eliminaras el componente 4 del diagrama (dejando solo 3), ¿aumentaría o disminuiría $E[Y]$?\n", "3. ¿Por qué el componente 5 aparece en serie con todo lo demás, y qué implica eso para la confiabilidad global del sistema (piensa en qué pasa si $X_5$ es muy pequeño)?\n", "4. Nuestra estimación de $\\bar Y(n)$ y $\\hat p_7$ no coincide exactamente con los valores del libro (4.33 y 0.19). ¿Es eso un error? Justifica usando lo que sabes sobre estimadores insesgados." ] }, { "cell_type": "markdown", "id": "df157d43", "metadata": {}, "source": [ "_Escribe aquí tus respuestas a las 4 preguntas anteriores._" ] }, { "cell_type": "markdown", "id": "4522fcfa", "metadata": {}, "source": [ "## 12. Conclusiones\n", "\n", "_Escribe 2-3 conclusiones sobre lo observado en este notebook._" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.11" } }, "nbformat": 4, "nbformat_minor": 5 }