{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2315e242",
   "metadata": {},
   "source": [
    "# CTMC con dos variables de estado — subsistemas de lavado (B) y llenado (C)\n",
    "## Extensión bivariada del Caso A (planta embotelladora)\n",
    "\n",
    "**Motivación.** En el ejercicio de la Clase 1 modelamos un solo subsistema de la planta (la Unidad de Lavado) con **una** variable de estado $X(t)\\in\\{0,1,2\\}$. Pero el paper de Gayathri (2025) en realidad describe **8 subsistemas** que evolucionan de forma simultánea — el modelo completo de 32 estados de ese paper se construye exactamente combinando las variables de estado de todos ellos. Aquí hacemos esa combinación de forma explícita, pero solo con **2 subsistemas** para que sea manejable: la Unidad de Lavado (B) y la Unidad de Llenado (C).\n",
    "\n",
    "Curso: Modelado de Sistemas bajo Incertidumbre · Universidad de los Andes · 2026-20\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24794d93",
   "metadata": {},
   "source": [
    "## Descripción del problema\n",
    "\n",
    "La línea de la planta embotelladora necesita que **ambos** subsistemas — lavado y llenado — estén operando para producir. Los dos subsistemas fallan y se reparan de forma **independiente** entre sí (una falla en uno no afecta al otro):\n",
    "\n",
    "- **Unidad de lavado (B)**: la misma del Caso A — una máquina principal y una de respaldo en frío. Tiene 3 posibles condiciones: Normal, Capacidad reducida, Falla total.\n",
    "- **Unidad de llenado (C)**: una sola máquina, sin respaldo. Tiene 2 posibles condiciones: Normal, Falla.\n",
    "\n",
    "**Dos variables de estado:**\n",
    "- $X_1(t)$ ....\n",
    "- $X_2(t)$: ...\n",
    "\n",
    "**Estado conjunto del sistema:** el par $(X_1(t),X_2(t))$. Como las dos unidades pueden estar en cualquier combinación de sus condiciones individuales, el **espacio de estados conjunto** es el producto cartesiano:\n",
    "\n",
    "$S = ... $\n",
    "\n",
    "es decir, **n estados combinados** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b794d84",
   "metadata": {},
   "outputs": [],
   "source": [
    "# importar librerías\n",
    "import numpy as np\n",
    "from jmarkov.ctmc import ctmc\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffe63aa3",
   "metadata": {},
   "source": [
    "## Parámetros del problema\n",
    "\n",
    "**Unidad de lavado (B)** — igual que en el Caso A:\n",
    "- $\\lambda_{1B}=0.02$/día, $\\lambda_{2B}=0.05$/día, $\\mu_{1B}=1$/día, $\\mu_{2B}=0.5$/día.\n",
    "\n",
    "**Unidad de llenado (C)** — una sola máquina, sin respaldo:\n",
    "- $\\lambda_C=0.03$ fallas/día (vida media 33.3 días).\n",
    "- $\\mu_C=0.8$ reparaciones/día (tiempo medio de reparación 1.25 días)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f7e839b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 1: define los parámetros de B y de C\n",
    "lda1_B, lda2_B, mu1_B, mu2_B = None, None, None, None\n",
    "lda_C, mu_C = None, None\n",
    "print(f'B: lda1={lda1_B}, lda2={lda2_B}, mu1={mu1_B}, mu2={mu2_B}')\n",
    "print(f'C: lda={lda_C}, mu={mu_C}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c408cfee",
   "metadata": {},
   "source": [
    "## Matrices generadoras individuales\n",
    "\n",
    "$$Q_B=\\begin{pmatrix}\\end{pmatrix} \\qquad Q_C=\\begin{pmatrix} \\end{pmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01df8ab0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 2: construye Q_B (3x3) y Q_C (2x2)\n",
    "Q_B = np.array([\n",
    "\n",
    "])\n",
    "Q_C = np.array([\n",
    "\n",
    "])\n",
    "print(f'Q_B=\\n{Q_B}')\n",
    "print(f'Q_C=\\n{Q_C}')\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert np.allclose(Q_B.sum(axis=1), 0) and np.allclose(Q_C.sum(axis=1), 0)\n",
    "print('OK: Q_B y Q_C son matrices generadoras válidas.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2a0384c",
   "metadata": {},
   "source": [
    "## Construir la matriz generadora conjunta — Forma 1: a mano\n",
    "\n",
    "Ordenamos los n estados combinados como $(b,c)$ con $b\\in $, $c\\in $, usando el índice $i=2b+c$. En cada estado, B y C pueden hacer una transición cada uno (nunca las dos a la vez, porque son procesos independientes y continuos):\n",
    "- si B cambia: la tasa es la misma que en $Q_B$ para esa fila de $b$, y $c$ no cambia.\n",
    "- si C cambia: la tasa es la misma que en $Q_C$ para esa fila de $c$, y $b$ no cambia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab2024a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 3: construye la matriz conjunta recorriendo los 6 estados (b,c)\n",
    "estados = [(b, c) for b in range(3) for c in range(2)]\n",
    "idx = {s: i for i, s in enumerate(estados)}\n",
    "n = len(estados)\n",
    "Q_manual = np.zeros((n, n))\n",
    "\n",
    "for (b, c) in estados:\n",
    "    i = idx[(b, c)]\n",
    "    # TODO: para cada b2 != b, si Q_B[b,b2] > 0, suma esa tasa en Q_manual[i, idx[(b2,c)]]\n",
    "    # TODO: para cada c2 != c, si Q_C[c,c2] > 0, suma esa tasa en Q_manual[i, idx[(b,c2)]]\n",
    "    # TODO: Q_manual[i,i] = -suma de la fila i (para que la fila sume 0)\n",
    "    pass\n",
    "\n",
    "print(f'estados={estados}')\n",
    "print(f'Q_manual=\\n{np.round(Q_manual,4)}')\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert np.allclose(Q_manual.sum(axis=1), 0)\n",
    "print('OK: Q_manual es una matriz generadora válida (6x6).')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c2c2122",
   "metadata": {},
   "source": [
    "## Construir la matriz generadora conjunta — Forma 2: suma de Kronecker\n",
    "\n",
    "Cuando dos CTMC evolucionan de forma **independiente**, su generador conjunto tiene una fórmula cerrada: la **suma de Kronecker**\n",
    "\n",
    "$$Q = Q_B\\otimes I_C + I_B\\otimes Q_C$$\n",
    "\n",
    "donde $I_B$ e $I_C$ son las matrices identidad del tamaño de $Q_B$ y $Q_C$, y $\\otimes$ es el producto de Kronecker (`np.kron` en Python). Esta fórmula da exactamente la misma matriz que construir el generador \"a mano\", pero en una sola línea — y es la manera estándar de combinar procesos independientes en cadenas de Markov."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "252d4058",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 4: construye Q usando la suma de Kronecker y compárala con Q_manual\n",
    "I_B = np.eye(3)\n",
    "I_C = np.eye(2)\n",
    "\n",
    "Q = None  # usa np.kron(Q_B, I_C) + np.kron(I_B, Q_C)\n",
    "print(f'Q (Kronecker)=\\n{np.round(Q,4)}')\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert np.allclose(Q, Q_manual), 'la suma de Kronecker debe coincidir con la construcción manual'\n",
    "print('OK: la suma de Kronecker coincide exactamente con Q_manual.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bcffccd",
   "metadata": {},
   "source": [
    "## Modelar con jmarkov\n",
    "\n",
    "Construye el modelo con `ctmc(Q)` y el vector inicial (arranca en $(0,0)$: ambas unidades en Normal)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4146db48",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 5: construye el modelo jmarkov y el vector inicial\n",
    "mc = ctmc(Q)\n",
    "alpha = np.zeros(6)\n",
    "alpha[idx[(0, 0)]] = None\n",
    "print(f'alpha={alpha}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92803202",
   "metadata": {},
   "source": [
    "## Disponibilidad conjunta\n",
    "\n",
    "La línea completa produce solo si **ambos** subsistemas están arriba: B no está en falla total ($b\\neq2$) **y** C no está en falla ($c\\neq1$). La disponibilidad conjunta es:\n",
    "\n",
    "$$Av_{conjunta}(t)=\\sum_{(b,c)\\,:\\,b\\neq2,\\ c\\neq1}\\pi_{(b,c)}(t)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "275f7c09",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 6: calcula pi(t) y la disponibilidad conjunta para varios tiempos\n",
    "idx_arriba = [idx[(b, c)] for b in [0, 1] for c in [0]]\n",
    "print(f'índices de estados \"ambas unidades arriba\": {idx_arriba}')\n",
    "\n",
    "tiempos = [1, 5, 20, 50]\n",
    "for t in tiempos:\n",
    "    pi_t = mc.transient_probabilities(t, alpha)\n",
    "    Av_t = None  # suma de pi_t en los índices \"arriba\"\n",
    "    print(f't={t}: Av_conjunta(t) = {Av_t:.5f}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f446d19c",
   "metadata": {},
   "source": [
    "## Distribución y disponibilidad estacionarias"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a60984f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 7: calcula la distribución estacionaria conjunta y la disponibilidad estacionaria\n",
    "ss = None  # usa mc.steady_state()\n",
    "print(f'pi (orden {estados}) = {np.round(ss,5)}')\n",
    "\n",
    "Av_ss = None\n",
    "print(f'Disponibilidad conjunta estacionaria = {Av_ss:.5f}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca023f84",
   "metadata": {},
   "source": [
    "## Conclusión\n",
    "\n",
    "_¿Por qué la disponibilidad conjunta es menor que la disponibilidad individual de cada subsistema? Verifica numéricamente si $Av_{conjunta}=Av_B\\times Av_C$ en estado estacionario (calcula $Av_B$ y $Av_C$ por separado, con sus propios `ctmc(Q_B)` y `ctmc(Q_C)`, y compara)._"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}