¿Se puede saber la resistencia del cemento sin esperar 28 días?

Qué está en juego

• Sin el modelo: diseño de mezcla lento, iterativo y costoso. Cada prueba física cuesta tiempo y material.
• Con el modelo: el ingeniero puede simular proporciones en un Excel y ver la resistencia estimada en segundos, priorizando qué mezclas vale la pena probar físicamente.
• Si el modelo falla: predicciones incorrectas pueden llevar a mezclas fuera de especificación — riesgo estructural y legal.

Para quién es este proyecto

Stack tecnológico

• Python 3.10+ — lenguaje principal
• Scikit-learn — Random Forest, SVR, KNN, PCA, VarianceThreshold, StandardScaler, GridSearchCV
• XGBoost — modelo de mayor ROI en benchmarking
• pandas / numpy / matplotlib — manipulación y visualización
• FastAPI + uvicorn — API REST para predicción
• ngrok — túnel para exponer la API desde Colab
• Google Sheets + Apps Script — interfaz de usuario para planta
• pickle — serialización del modelo entrenado

¿Qué está hecho el cemento?

Para trabajar con datos de cemento no necesitas ser químico, pero sí necesitas entender el dominio mínimo. El cemento Tipo ICo es una mezcla de componentes con proporciones variables. Los principales son clínker conforme, clínker no conforme, yeso, piedra caliza, arcilla calcinada, escoria y varios óxidos químicos (SiO₂, Al₂O₃, Fe₂O₃, etc.).

La hipótesis central del proyecto — definida junto con el equipo de I+D de la empresa — era que el yeso tiene una influencia directa y medible en la resistencia a la compresión. No cualquier componente: específicamente el yeso. Eso guió desde qué variables mirar primero hasta cómo interpretar los resultados del modelo.

Definición formal del problema

• Variable objetivo (y): Resistencia a la compresión en PSI a 28 días de fraguado.
• Variables de entrada (X): Composición química del cemento — proporciones en porcentaje de cada componente.
• Tipo de tarea: Regresión supervisada.
• Métricas válidas: RMSE (penaliza errores grandes, relevante porque subestimar resistencia es peligroso), MAE (error promedio absoluto en PSI), MAPE (error relativo en %).
• Métricas que NO aplican: accuracy, precision, recall, F1 — son métricas de clasificación.
• Nivel de error aceptable: En la literatura industrial, un MAPE de 5–10% sobre valores de resistencia es considerado práctico para aplicaciones de control de calidad. Valores mayores al 15% probablemente no justifican reemplazar las pruebas físicas.

¿Qué predicciones son peligrosas?

No todos los errores son iguales. En este problema, subestimar la resistencia (predecir 3,800 PSI cuando la real es 4,500 PSI) lleva a descartar una mezcla buena — desperdicio de material y tiempo. Pero sobreestimar la resistencia (predecir 4,500 PSI cuando la real es 3,800 PSI) es peor: podría llevar a aprobar una mezcla que no cumple especificaciones, con consecuencias estructurales y legales.

Esto tiene implicaciones en cómo interpretas las métricas: un modelo con RMSE bajo pero sesgado hacia sobreestimar es más peligroso que uno con mayor error pero sin sesgo sistemático.

Sobre los datos públicos: cuándo puedes usarlos como proxy

Los datos reales de la empresa son confidenciales. Para replicar el proceso, puedes usar dos datasets públicos como proxy:

Análisis cualitativo antes del cuantitativo

El primer paso no es correr código. Es sentarte con alguien que entiende el dominio — en este caso el equipo de Ingeniería Química e I+D — y clasificar cada variable antes de ver los números:

• Verde — variable altamente útil, se espera que tenga relación directa con la resistencia.
• Ámbar — incierta, puede o no ser relevante, requiere verificación cuantitativa.
• Rojo — irrelevante a priori, o derivada de otras variables (riesgo de leakage).

Este paso te ahorra entrenar modelos con variables que ya sabes que son ruido. También te da una hipótesis que el modelo deberá confirmar o refutar.

Carga y limpieza inicial

El dataset es un Excel con 31 columnas de composición química y 4 columnas target. La columna Clínker I Total (%) se elimina porque es combinación lineal de otras columnas — incluirla sería introducir redundancia perfecta. Arcilla Calcinada (%) también se descarta porque todos sus valores son cero: varianza nula, sin información.

import pandas as pd

data = pd.read_excel("data/composito_cemento_tipo_ico.xlsx", header=1)
data = data.dropna()

# Atributos: columnas 2 a 33, excluyendo columnas redundantes
df_atributes = data.iloc[:, 2:34].copy()
del df_atributes["Clínker I Total (%)"]        # combinación lineal → leakage
df_atributes.drop(["Arcilla Calcinada (%)"], axis=1, inplace=True)  # varianza 0

# Target: resistencia a 28 días
y = data.iloc[:, -5]

print(f"Features: {df_atributes.shape[1]}")   # debe ser ~29
print(f"Muestras: {len(y)}")
print(f"\nDistribución del target (PSI):")
print(y.describe())

Análisis cuantitativo de variables: VarianceThreshold y correlación

from sklearn.feature_selection import VarianceThreshold
import pandas as pd

# Eliminar features con varianza muy baja (casi constantes)
sel = VarianceThreshold(threshold=(.8 * (1 - .8)))
df_filtrado = sel.fit_transform(df_atributes)

# Ver qué variables sobrevivieron
features_selected = df_atributes.columns[sel.get_support()]
features_dropped = df_atributes.columns[~sel.get_support()]

print(f"Features originales: {df_atributes.shape[1]}")
print(f"Features tras VarianceThreshold: {df_filtrado.shape[1]}")
print(f"\nEliminadas por baja varianza: {list(features_dropped)}")

import matplotlib.pyplot as plt

corr_target = df_atributes.corrwith(y).sort_values(ascending=False)

plt.figure(figsize=(10, 6))
corr_target.plot.barh()
plt.title("Correlación de Pearson de cada feature con RC 28 días (PSI)")
plt.xlabel("Correlación")
plt.axvline(x=0, color="gray", linestyle="--", linewidth=0.8)
plt.tight_layout()
plt.show()

print("Top 5 correlacionadas positivamente:")
print(corr_target.head())
print("\nTop 5 correlacionadas negativamente:")
print(corr_target.tail())

Reducción de dimensiones con PCA

Con ~29 variables altamente correlacionadas entre sí, PCA reduce la dimensionalidad conservando la varianza explicada. Se entrena PCA completo primero para ver cuántos componentes son necesarios, luego se aplica con ese número.

import numpy as np
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

pca_full = PCA()
pca_full.fit(df_atributes)

plt.figure(figsize=(10, 5))
plt.plot(np.cumsum(pca_full.explained_variance_ratio_), marker="o", markersize=4)
plt.xlabel("Número de componentes")
plt.ylabel("Varianza acumulada explicada")
plt.axhline(y=0.95, color="red", linestyle="--", label="95% varianza")
plt.legend()
plt.title("PCA — Varianza acumulada")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Aplicar PCA con el número óptimo (ej. 10 componentes → ~95% varianza)
pca = PCA(n_components=10)
df_reducido = pca.fit_transform(df_atributes)
print(f"Varianza explicada con 10 componentes: {pca.explained_variance_ratio_.sum():.2%}")

Split y normalización

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

X = df_atributes.copy()

X_train, X_temp, y_train, y_temp = train_test_split(X, y, test_size=0.2, random_state=42)
X_val, X_test, y_val, y_test = train_test_split(X_temp, y_temp, test_size=0.5, random_state=42)

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_val_scaled   = scaler.transform(X_val)
X_test_scaled  = scaler.transform(X_test)

print(f"Train: {X_train_scaled.shape} | Val: {X_val_scaled.shape} | Test: {X_test_scaled.shape}")

¿Por qué estos cuatro modelos?

• Regresión Lineal — baseline obligatorio. Si los demás modelos no lo superan claramente, hay un problema en los datos o en el pipeline.
• Random Forest — robusto a correlaciones entre features, provee importancia de características de forma nativa. Buena interpretabilidad para el ingeniero.
• XGBoost — mejor precisión en benchmarking (menor RMSE), pero ligeramente más complejo de explicar. Máximo ROI relativo.
• Red Neuronal (MLP) — incluida para comparación. En este tipo de dato estructurado y con pocos registros, raramente supera a los ensembles.

¿Por qué no SVR? Se incluyó en iteraciones tempranas del proyecto pero su sensibilidad a los hiperparámetros C y epsilon lo hace costoso de tunear sin ganancia clara frente a Random Forest. Se reemplazó por XGBoost en la versión final.

Sintonización del modelo seleccionado

from sklearn.model_selection import GridSearchCV
from xgboost import XGBRegressor

param_grid = {
    "n_estimators":     [100, 200, 300],
    "max_depth":        [3, 5, 7],
    "learning_rate":    [0.05, 0.1, 0.2],
    "subsample":        [0.8, 1.0],
    "colsample_bytree": [0.8, 1.0],
}

grid_search = GridSearchCV(
    XGBRegressor(random_state=42, verbosity=0),
    param_grid,
    cv=5,
    scoring="neg_root_mean_squared_error",
    n_jobs=-1,
    verbose=1
)

grid_search.fit(X_train_scaled, y_train)

print(f"Mejores parámetros: {grid_search.best_params_}")
print(f"Mejor RMSE CV: {-grid_search.best_score_:.2f} PSI")

Guardar modelo y scaler

import pickle

mejor_modelo = grid_search.best_estimator_
mejor_modelo.fit(X_train_scaled, y_train)

with open("models/modelo_cemento_xgb.pkl", "wb") as f:
    pickle.dump(mejor_modelo, f)

with open("models/scaler_cemento.pkl", "wb") as f:
    pickle.dump(scaler, f)

# También guarda los nombres de columnas: el orden importa en producción
import json
with open("models/columnas_modelo.json", "w") as f:
    json.dump(list(df_atributes.columns), f)

print("Modelo, scaler y columnas guardados.")

Evaluación en el conjunto de test

from sklearn.metrics import mean_squared_error, mean_absolute_error
import numpy as np
import matplotlib.pyplot as plt

y_pred = mejor_modelo.predict(X_test_scaled)
y_true = y_test.values

mae  = mean_absolute_error(y_true, y_pred)
rmse = np.sqrt(mean_squared_error(y_true, y_pred))
mape = np.mean(np.abs((y_true - y_pred) / y_true)) * 100

print(f"MAE:  {mae:.2f} PSI  — error promedio en PSI")
print(f"RMSE: {rmse:.2f} PSI  — penaliza errores grandes")
print(f"MAPE: {mape:.2f}%    — error relativo promedio")

# Gráfica: valor real vs predicho
plt.figure(figsize=(8, 6))
plt.scatter(y_true, y_pred, alpha=0.6, edgecolors="none")
plt.plot([y_true.min(), y_true.max()], [y_true.min(), y_true.max()],
         "r--", linewidth=1.5, label="Predicción perfecta")
plt.xlabel("Resistencia real (PSI)")
plt.ylabel("Resistencia predicha (PSI)")
plt.title("Real vs. Predicho — XGBoost sintonizado")
plt.legend()
plt.tight_layout()
plt.show()

Comparativa final de modelos

El ROI relativo se estima con la fórmula: ROI = [(RMSE_base − RMSE_modelo) × C_error + Ahorro_tiempo + Ahorro_operativo − Costo_modelo] / Costo_modelo. El factor C_error es el costo de cada unidad de error en PSI para el negocio — definido con el cliente.

Importancia de características — ¿el yeso realmente importa?

import pandas as pd
import matplotlib.pyplot as plt

importancias = pd.Series(
    mejor_modelo.feature_importances_,
    index=df_atributes.columns
).sort_values(ascending=False)

plt.figure(figsize=(10, 6))
importancias.head(15).plot.bar(color="#00b76c")
plt.title("Top 15 features por importancia — XGBoost")
plt.ylabel("Importancia relativa")
plt.xticks(rotation=45, ha="right")
plt.tight_layout()
plt.show()

print("Top 5 features:")
print(importancias.head())

¿Cuándo NO confiar en el modelo?

• Cuando cambia la composición fuera del rango entrenado: si introduces proporciones de yeso o clínker que no existían en los datos de entrenamiento, las predicciones son extrapolaciones no validadas.
• Cuando cambian condiciones no modeladas: temperatura de cocción, calidad del agua usada en planta, cambio de proveedor de materia prima. El modelo asume que la composición química es suficiente para explicar la resistencia — si esas condiciones externas cambian, el modelo pierde validez.
• Cuando el input es incorrecto: si el ingeniero ingresa proporciones que no suman correctamente o hay errores de medición, la predicción será inválida. El modelo no tiene mecanismo para detectarlo.

Arquitectura del prototipo

El flujo completo tiene cuatro componentes:

• Google Sheets UI: el usuario ingresa los datos de composición y ve la predicción en la misma celda.
• Apps Script Function (predecirRC()): convierte los datos de la hoja a JSON, hace el request HTTPS al servidor y escribe la predicción de vuelta.
• Túnel ngrok: mapea https://xxx.ngrok.dev al puerto local 8000 donde corre FastAPI.
• Servidor (Google Colab VM): carga el modelo .pkl y el scaler, recibe el JSON, normaliza, predice y retorna el resultado.

Cómo correr el prototipo completo

# En Google Colab, ejecutar en orden:

# 1. Instalar dependencias
!pip install fastapi uvicorn pyngrok xgboost -q

# 2. Escribir el archivo api/main.py (ver código arriba)
# ... (usa %%writefile api/main.py)

# 3. Levantar el servidor en background
import subprocess
proc = subprocess.Popen(
    ["uvicorn", "api.main:app", "--host", "0.0.0.0", "--port", "8000"],
    stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL
)

# 4. Abrir túnel ngrok
from pyngrok import ngrok
ngrok.set_auth_token("TU_TOKEN_AQUI")  # registrarse gratis en ngrok.com
tunnel = ngrok.connect(8000)
print(f"API disponible en: {tunnel.public_url}")
# 👉 Copia esta URL y pégala en tu Apps Script