Simulación resuelta: del balance a una conclusión comprobable

Justificar una conclusión con una simulación

Los capítulos 21 a 29 separaron problemas que suelen aparecer mezclados cuando se construye una simulación, como la representación del estado, la organización del código y la visualización, la integración numérica y la rigidez, la precisión finita, la exploración paramétrica y la verificación. Este capítulo los reúne en una sola cadena de trabajo. Sin añadir métodos numéricos nuevos, sigue un modelo desde la formulación de una pregunta hasta una conclusión cuyo alcance pueda auditarse.

El ejemplo es deliberadamente sencillo. Una carga acumulada recibe entradas a una tasa constante y pierde unidades a una tasa proporcional a su tamaño. Podría tratarse de un inventario o del material de un compartimento. Usaremos la interpretación de una carga abstracta de trabajo para mantener la conexión con los sistemas organizacionales, aunque el modelo no pretende describir un servicio clínico real. Como la ecuación tiene solución analítica, podremos distinguir los errores numéricos de las limitaciones estructurales del propio modelo.

La pregunta es concreta. Si la carga parte de 120 unidades, ingresan ocho unidades por día y la salida instantánea equivale al diez por ciento de la carga por día, ¿cuánto tarda el sistema en bajar hasta 90? Después examinaremos cuánto depende la respuesta del parámetro de salida y qué ocurre si la salida proporcional se reemplaza por una capacidad fija.

La secuencia obliga a mantener separadas la definición del problema y el modelo matemático que la formaliza, así como la solución numérica de ese modelo y la interpretación del resultado. Una simulación científicamente útil requiere coherencia entre esas capas. Una curva visualmente convincente no sustituye a ninguna de ellas.

Delimitación del sistema y ecuación de balance

Sea $W(t)$ la cantidad acumulada de trabajo en el instante $t$, medida en unidades de carga. Definimos una entrada constante $\lambda=8$ unidades por día y una salida proporcional a la carga, $O(t)=kW(t)$, con $k=0{,}10\ \text{día}^{-1}$. La ecuación de balance es entonces

Las unidades dan una primera comprobación. El término $dW/dt$ tiene unidades de carga por día, igual que $\lambda$. Como $k$ tiene unidades de día inverso, $kW$ vuelve a tener unidades de carga por día, de modo que la ecuación es dimensionalmente coherente.

La frontera del modelo es igualmente importante. El modelo no incluye tipos distintos de trabajo ni prioridades. Tampoco incluye una capacidad que dependa del personal disponible, llegadas aleatorias o decisiones adaptativas. No distingue individuos, porque $W$ es una variable agregada. Si alguno de esos mecanismos es relevante para la pregunta empírica, la ecuación puede ser insuficiente aunque se resuelva con exactitud.

Antes de usar un computador pueden derivarse varias consecuencias. El equilibrio satisface $0=\lambda-kW_{\mathrm{eq}}$, de modo que

Como la condición inicial es mayor que 80, cabe esperar una disminución monótona hacia ese equilibrio. El umbral operativo que queremos estudiar es $H=90$. Dado que $W_{\mathrm{eq}}<H<W(0)$, el modelo predice que el umbral se cruzará en un tiempo finito. Estas expectativas, obtenidas sin simular, servirán después para verificar el cálculo.

Solución analítica y tiempo exacto hasta el evento

La ecuación es lineal y puede resolverse directamente. Con condición inicial $W_0$, la solución es

que con nuestros parámetros da $W(t)=80+40e^{-0{,}1t}$. El tiempo exacto en que $W=90$ se obtiene de $90=80+40e^{-0{,}1t}$, es decir, de $e^{-0{,}1t}=1/4$, lo que da

Este resultado cumple dos funciones pedagógicas. Responde exactamente la pregunta dentro del modelo y da un estándar externo para comprobar la implementación numérica. Si el programa entregara un cruce a los 18 días, habría un problema en la integración, en la detección del evento o en la especificación de los parámetros. Esa discrepancia no debería interpretarse como una propiedad del sistema.

La solución analítica permite también comprobar límites. En $t=0$ da 120 y, cuando $t\to\infty$, el término exponencial desaparece y $W\to80$. La derivada inicial es

y la magnitud de la disminución se reduce a medida que el sistema se acerca al equilibrio, de modo que la trayectoria no es una recta.

Especificación computacional y criterio de evento

Aunque la solución exacta se conoce, resolveremos el problema numéricamente, porque el flujo de trabajo debe poder trasladarse a modelos que no la tengan. Usaremos un integrador adaptativo para problemas de valor inicial y definiremos explícitamente la función de evento.

from scipy.integrate import solve_ivp

lam = 8.0
k = 0.10
W0 = 120.0
H = 90.0

def rhs(t, y):
    return [lam - k*y[0]]

def evento(t, y):
    return y[0] - H

evento.terminal = True
evento.direction = -1

sol = solve_ivp(
    rhs,
    (0.0, 50.0),
    [W0],
    rtol=1e-10,
    atol=1e-12,
    events=evento,
    dense_output=True,
)

La función solve_ivp integra numéricamente sistemas de ecuaciones diferenciales ordinarias con condiciones iniciales y permite localizar eventos definidos como ceros de funciones continuas (SciPy Developers, s. f.). El instante del evento se estima sobre la solución interpolada entre dos pasos (Hairer et al., 1993). La elección de las tolerancias pertenece a la especificación numérica y es independiente de la ecuación científica.

El evento se define como $g(W)=W-H$. El atributo direction, fijado en menos uno, indica que interesa el cruce desde valores positivos hacia negativos, coherente con la pregunta. La opción terminal = True detiene la integración cuando se identifica ese cruce. Estas decisiones parecen menores, pero forman parte del significado del resultado. En otro modelo, un detector sin dirección podría registrar un cruce de regreso. Un evento mal definido, por su parte, podría localizar una cantidad distinta de la pretendida.

Con estas condiciones y SciPy 1.17.1, la simulación obtiene $t_H^{\text{num}}=13{,}8629436123$ días, frente al valor exacto $t_H=13{,}8629436112$ días. La diferencia es de unos $1{,}1\times10^{-9}$ días, irrelevante para la pregunta planteada. Esa magnitud pequeña no demuestra, sin embargo, que el modelo sea empíricamente correcto. Solo muestra que, en este caso, la solución numérica reproduce con gran precisión la solución matemática especificada (figura 30.1).

Verificación mínima antes de interpretar

Una conclusión auditada requiere algo más que comparar una cifra final. En este modelo pueden reunirse varias verificaciones pequeñas. La consistencia dimensional ya se comprobó, porque todos los términos del balance tienen unidades de carga por día. El cálculo debe devolver exactamente la condición inicial $W(0)=120$. Si se parte de $W=80$, la derivada debe anularse dentro de la precisión numérica.

La monotonía da otra prueba. Con el estado inicial usado, $W(t)>80$ y, por tanto, $dW/dt<0$, de modo que una trayectoria que aumentara antes de llegar al equilibrio contradiría el modelo. La solución exacta permanece, además, positiva. Un método que produjera valores negativos en este escenario obligaría a revisar el paso o la implementación. La trayectoria numérica y el tiempo del evento pueden compararse, por último, con sus expresiones analíticas.

La convergencia completa el conjunto. Con Euler explícito, el error en un tiempo fijo debe disminuir al reducir el paso. El cuadro 30.1 muestra que en $t=10$ el error se divide aproximadamente por dos cada vez que el paso se divide por dos, como corresponde a un método de primer orden. Euler no es el integrador recomendado para este problema, pero su comportamiento esperado permite comprobar la implementación.

Cuadro 30.1. Método de Euler en $t=10$ para el modelo base, cuyo valor exacto es $W(10)=80+40e^{-1}\approx94{,}715178$. El cociente divide el error de la fila anterior por el de la fila actual.

Estas pruebas no son equivalentes. Una implementación puede conservar el equilibrio y calcular mal los transitorios. Otra puede coincidir con un valor final por compensación de errores. La fuerza de la verificación procede de combinar propiedades parcialmente independientes, como vimos en el capítulo 29.

Dependencia de la conclusión respecto del parámetro

La tasa $k$ controla tanto el equilibrio como la velocidad de respuesta. Si $\lambda=8$ se mantiene fijo, $W_{\mathrm{eq}}=8/k$. La pregunta sobre el tiempo hasta 90 tiene respuesta finita solo cuando el equilibrio está por debajo de 90, es decir, cuando $8/k<90$, lo que exige

Este valor es un umbral estructural del modelo respecto de nuestra definición de evento. Con $k=0{,}08$ el equilibrio es 100, de modo que la carga baja desde 120 sin alcanzar nunca 90. Con $k=0{,}09$ el equilibrio es de unos 88,89 y el cruce existe, aunque tarda unos 37,02 días. Con $k=0{,}15$ ocurre en unos 3,99 días (cuadro 30.2).

Cuadro 30.2. Equilibrio y tiempo hasta $H=90$ para distintos valores de $k$, con $\lambda=8$ y $W(0)=120$.

La figura 30.2 muestra que modificar $k$ cambia la pendiente inicial y, además, el estado al que tiende el sistema.

La expresión exacta del tiempo de cruce es

válida cuando $\lambda/k<H<W_0$. Cuando $k$ se acerca desde arriba al valor crítico $\lambda/H$, el numerador dentro del logaritmo tiende a cero y el tiempo crece sin cota (figura 30.3). Una extrapolación lineal de tiempos observados lejos del umbral describe mal esa transición.

La sensibilidad puede examinarse también en el equilibrio. Como $W_{\mathrm{eq}}=\lambda/k$, se tiene

Estas derivadas tienen unidades, de modo que su magnitud depende de las escalas elegidas. Las elasticidades adimensionales son

En este modelo, un aumento relativo pequeño de la entrada produce el mismo cambio relativo en el equilibrio, mientras que un aumento relativo pequeño de $k$ lo reduce en la misma proporción. Esa afirmación local no se traslada al tiempo del evento. La elasticidad de $t_H$ respecto de $k$ vale unos $-5{,}3$ en $k=0{,}10$ y unos $-24$ en $k=0{,}09$. Cerca de $k=\lambda/H$, una pequeña modificación de $k$ puede cambiar el tiempo de cruce en decenas de días o eliminar el cruce finito. La sensibilidad depende de cuál salida se observa.

El principio metodológico reaparecerá a menudo. La sensibilidad de un modelo depende de la pregunta, porque un parámetro puede afectar modestamente una variable y de manera extrema otra, sobre todo cuando esta última se define mediante umbrales o tiempos de primera llegada. Todo análisis de sensibilidad debe especificar, por eso, el parámetro, la salida, el dominio y la escala de la perturbación.

También hay que distinguir la sensibilidad de la incertidumbre. Que el resultado sea muy sensible a $k$ nada dice de cuánto se desconoce $k$. Si el parámetro se midiera con enorme precisión, la incertidumbre inducida podría seguir siendo pequeña. A la inversa, un parámetro de efecto moderado pero muy mal determinado puede contribuir mucho a la incertidumbre final. En las aplicaciones reales hay que combinar ambas preguntas, es decir, cuánto cambia el resultado cuando cambia el parámetro y qué valores del parámetro son compatibles con la evidencia disponible. Esa combinación pertenece al análisis de incertidumbre, que estudiaremos más adelante con herramientas formales.

Barrido, calibración e identificación

La exploración anterior responde a una pregunta condicional sobre las consecuencias del modelo si $k$ toma determinados valores. Nada se estimó a partir de datos. Tampoco se estableció que $k$ sea constante ni que una salida real sea proporcional a la carga.

La distinción es central. Un barrido paramétrico explora la estructura de un modelo, mientras que la calibración busca relacionar los parámetros con observaciones. La identificación pregunta si los datos disponibles permiten determinar esos parámetros de manera suficientemente única. Las tres operaciones pueden usar el mismo símbolo $k$, pero tienen compromisos epistémicos distintos.

Supongamos que una organización observa la disminución de la carga durante diez días. Podría ajustar un valor de $k$ que reprodujera esa serie. Otras estructuras, como una capacidad fija o una salida saturable, podrían producir trayectorias parecidas en un intervalo limitado. Si la serie de diez días procediera de una capacidad fija de 12 unidades por día, el ajuste por mínimos cuadrados del modelo proporcional daría $k\approx0{,}116$ con un error cuadrático medio de 2,5 unidades y predeciría el cruce de 90 a los 7,6 días, muy cerca de los 7,5 correctos. Esos dos modelos discrepan, sin embargo, por completo a largo plazo, porque uno tiende a un equilibrio de unas 69 unidades y el otro vacía la carga en 30 días. Un buen ajuste local no demuestra que la ley proporcional sea el mecanismo correcto (Sterman, 2000). El problema reaparecerá cuando estudiemos la identificabilidad y la comparación de modelos.

Incertidumbre numérica e incertidumbre estructural

Hasta aquí el error numérico es diminuto. Podría resultar tentador informar el tiempo de 13,86294 días con muchas cifras y presentar esa precisión como una propiedad del problema, lo que sería un error conceptual. Para verlo, sustituyamos la salida proporcional por una capacidad fija de 12 unidades por día mientras exista carga. El balance pasa a ser $dW/dt=8-12=-4$ mientras $W>0$. Desde 120, el umbral de 90 se alcanza en $(120-90)/4=7{,}5$ días.

El modelo proporcional predice 13,86 días y el de capacidad fija, 7,5 (figura 30.4). La diferencia entre modelos es de varios días, mientras que el error del integrador para el modelo proporcional es del orden de una milmillonésima de día.

De esto no se sigue que la capacidad fija sea mejor, porque ambos son modelos estilizados. La comparación muestra una jerarquía de incertidumbres. Reducir el error numérico más allá de cierto punto no compensa ignorar una incertidumbre estructural dominante.

Por eso hay que documentar la forma del modelo junto con sus parámetros y tolerancias. La norma NASA-STD-7009B y su manual de implementación tratan la credibilidad de los modelos y simulaciones como un problema que incluye verificación, validación, procedencia de las entradas, incertidumbre, robustez y gestión del proceso (National Aeronautics and Space Administration, 2024, 2026). Según esos documentos, la credibilidad no puede deducirse solo de la resolución numérica.

Qué significa validar este modelo

La verificación pregunta si se resuelven correctamente las ecuaciones elegidas. La validación exige confrontar las consecuencias relevantes del modelo con el sistema de interés, para el uso previsto (Oberkampf y Roy, 2010; Sargent, 2013). En nuestro ejemplo no bastaría con observar que una curva real desciende.

Habría que examinar, entre otras cosas, si la entrada puede aproximarse por una tasa constante durante el período de interés y si la salida aumenta aproximadamente con la carga. También habría que averiguar si $k$ permanece estable, si hay retardos y si la variable agregada $W$ conserva la información necesaria para la decisión. La validación depende del propósito, porque un modelo que aproxima bien el promedio semanal podría ser insuficiente para predecir máximos diarios o tiempos individuales.

Tampoco existe una validación absoluta de una representación compleja. La pregunta práctica es si la evidencia disponible respalda su uso dentro de un dominio especificado y con una incertidumbre compatible con la decisión. Esa formulación evita dos extremos, el de tratar todo modelo como una copia verdadera del sistema y el de concluir que ninguna evaluación empírica es posible porque todo modelo simplifica.

Un informe reproducible del resultado

Una simulación auditable debe conservar información suficiente para reconstruir qué se calculó. En este ejemplo, el registro mínimo incluiría la ecuación del modelo con sus unidades y la definición de la variable de estado, la condición inicial $W_0=120$, los parámetros $\lambda=8$ y $k=0{,}10$ y el umbral $H=90$. Añadiría el intervalo temporal, el integrador y sus tolerancias, la definición y la dirección del evento y la versión del software. Lo completarían las pruebas de verificación realizadas y el código necesario para reproducir las figuras y los cuadros.

La ejecución de este capítulo usó SciPy 1.17.1. El integrador hizo 182 evaluaciones de la función en el escenario base con las tolerancias especificadas. Ese número caracteriza esta ejecución y puede cambiar con la versión, el método o las tolerancias, así que no debe confundirse con una propiedad del fenómeno. Registrar metadatos de este tipo ayuda a reconstruir el cálculo sin mezclar los detalles de implementación con los parámetros científicos.

El archivo complementario cap30_simulacion.py genera las cuatro figuras, los cuadros de escenarios y de convergencia y un resumen de verificación. Así los resultados del texto no dependen de copiar cifras a mano desde una sesión interactiva.

La conclusión auditada

Ahora puede formularse una conclusión proporcionada a la evidencia disponible.

Dentro del modelo $dW/dt=8-0{,}1W$, con $W(0)=120$, la carga disminuye monótonamente hacia un equilibrio de 80 unidades y cruza el umbral de 90 a los 13,86294 días. La solución numérica reproduce el tiempo analítico con un error aproximado de $1{,}1\times10^{-9}$ días bajo las tolerancias especificadas. El resultado depende fuertemente del parámetro $k$ y de la forma funcional supuesta para la salida, por lo que no debe interpretarse como predicción de un sistema real sin estimación de los parámetros y validación de la estructura del modelo.

La conclusión separa la consecuencia matemática de la adecuación numérica que la respalda y termina limitando la inferencia empírica. Esas partes son necesarias para evitar que la precisión computacional se convierta en una falsa apariencia de certeza científica.

Una conclusión menos cuidadosa, como «el sistema tardará 13,86 días en bajar de 90», elimina precisamente la información que permite juzgar su alcance. Los condicionantes de la conclusión especifican el conjunto de situaciones en que pretende ser válida.

Auditoría final del proceso

El flujo de trabajo puede resumirse mediante las preguntas del cuadro 30.3.

Cuadro 30.3. Auditoría final del proceso de simulación.

La auditoría no garantiza que el modelo sea adecuado. Su función es hacer visibles los puntos en que una afirmación puede fallar y separar los ya comprobados de los que todavía requieren evidencia externa.

Consolidación de la tercera parte

Esta parte del tomo comenzó en el capítulo 21 con una regla discreta pequeña y termina con una cadena completa de razonamiento computacional. El modelo especifica relaciones matemáticas y el código implementa operaciones destinadas a estudiarlas, de modo que verificar el código no valida el modelo. Una trayectoria simulada pertenece al modelo, mientras que los datos pertenecen a un procedimiento de observación. Compararlos exige definir el observable que los conecta.

Una figura puede revelar estructura o esconderla, porque las escalas, el muestreo, la normalización y los promedios forman parte del análisis. El paso temporal y el redondeo afectan la solución calculada. Los parámetros inciertos y la forma funcional afectan, en cambio, la relación entre el modelo y el fenómeno. Ambas clases de error no deben reunirse en una categoría vaga de «incertidumbre».

Explorar valores hipotéticos de un parámetro no equivale a estimarlo a partir de observaciones. Un resultado que no ocurre dentro de una ventana simulada puede ocurrir después, porque el horizonte forma parte de la pregunta. La verificación pregunta si el modelo especificado se resuelve con suficiente fidelidad, mientras que la validación pregunta si el modelo es adecuado para un uso empírico determinado. Añadir decimales a una salida, por último, no aumenta la credibilidad científica, porque una cifra puede calcularse con enorme precisión dentro de una representación equivocada.

Con estas herramientas, el computador deja de ser una caja que recibe ecuaciones y devuelve gráficos. Pasa a ser un instrumento cuyo comportamiento también debe verificarse y documentarse. Esa disciplina será necesaria en la cuarta parte, que a partir del capítulo 31 estudia las ecuaciones diferenciales, la estabilidad y la geometría del plano.

Ejercicios

Ejercicio 30.1. (a) Verifique que $W(t)=80+40e^{-0{,}1t}$ satisface la ecuación de balance y la condición inicial. (b) Calcule $t_H$ para $H=100$ y para $H=85$. (c) ¿Qué ocurre si se pide el tiempo hasta $H=80$?

Ejercicio 30.2. (a) Calcule el orden observado de Euler con los errores del cuadro 30.1 para $h=0{,}25$ y $h=0{,}125$. (b) Muestre que la aproximación de Euler en $t=10$ es $W_N=80+40(1-0{,}1h)^{10/h}$ y compruebe el valor para $h=2$. (c) ¿Con qué pasos produciría Euler aproximaciones que alternan alrededor del equilibrio y con cuáles divergiría?

Ejercicio 30.3. (a) Calcule numéricamente $dt_H/dk$ en $k=0{,}10$ y en $k=0{,}09$ mediante diferencias centradas. (b) Calcule las elasticidades correspondientes y compárelas con la elasticidad del equilibrio. (c) Explique por qué la sensibilidad depende de la salida elegida.

Ejercicio 30.4. Una organización registra durante diez días una carga que baja linealmente de 120 a 80, como produciría una capacidad fija de 12 unidades por día con una entrada de 8. (a) Ajuste por mínimos cuadrados el valor de $k$ del modelo proporcional con $\lambda=8$ y $W(0)=120$. (b) Compare el tiempo hasta 90 que predice cada modelo. (c) Compare sus predicciones a largo plazo.

Ejercicio 30.5. Redacte la conclusión auditada del capítulo para el modelo de capacidad fija, con la misma estructura de tres partes.

Ejercicio 30.6. Con el archivo complementario del capítulo, (a) cambie las tolerancias a rtol=1e-6 y atol=1e-9 y compare el tiempo del evento y el número de evaluaciones. (b) ¿Cambia la conclusión auditada? (c) ¿Qué información del registro permitiría a otra persona detectar ese cambio?

Referencias

Hairer, E., Nørsett, S. P. y Wanner, G. (1993). Solving ordinary differential equations I: Nonstiff problems (2.ª ed. rev.). Springer.

National Aeronautics and Space Administration. (2024). Standard for models and simulations (NASA-STD-7009B). NASA Office of the Chief Engineer. https://standards.nasa.gov/sites/default/files/standards/NASA/B/1/NASA-STD-7009B-Final-3-5-2024.pdf

National Aeronautics and Space Administration. (2026). NASA handbook for models and simulations: An implementation guide for NASA-STD-7009B (NASA-HDBK-7009B). NASA Office of the Chief Engineer. https://standards.nasa.gov/standard/nasa/nasa-hdbk-7009

Oberkampf, W. L. y Roy, C. J. (2010). Verification and validation in scientific computing. Cambridge University Press.

Sargent, R. G. (2013). Verification and validation of simulation models. Journal of Simulation, 7(1), 12–24. https://doi.org/10.1057/jos.2012.20

SciPy Developers. (s. f.). scipy.integrate.solve_ivp. Recuperado el 19 de septiembre de 2026, de https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Sterman, J. D. (2000). Business dynamics: Systems thinking and modeling for a complex world. Irwin/McGraw-Hill.