Evaluación de riesgo en el cronograma con VBA Excel

Este tema es más aplicable a la gestión de proyectos, específicamente para la gestión del riesgo en el plazo donde uno de los elementos a controlar es la certeza del cronograma que se elabora.

Dependiendo del nivel de descomposición del trabajo para cada actividad se estima un tiempo de duración y, de la culminación de cada actividad (o de parte de ella), depende el inicio de una o más actividades(según el secuenciado de actividades).

La definición de la duración de las actividades se hace mediante Estimación Análoga, el Juicio de Expertos, Estimación por 3 valores, etc. esto último se enmarca dentro de la técnica conocida como PERT-CPM. Para el conjunto de actividades y su correspondiente secuenciado, es posible obtener estimaciones en base a los cálculos de las varianzas de la ruta crítica que nos sirvan para una evaluación del cronograma. En este caso se asume que la duración del proyecto se ajusta a una distribución normal.

El método PERT utiliza tres estimaciones para definir un rango aproximado de duración de una actividad.
  • duración más probable (tm)
  • duración optimista (to)
  • duración pesimista (tp)
La duración esperada se calcula con:  te = (to + 4 tm + tp) / 6

Esto es aplicable a las actividades para las que se espera que su duración tenga un comportamiento del tipo PERT. Es posible que algunos procesos no tienen un comportamiento PERT sino que podrían tener un comportamiento del tipo triangular para el cual:  te = (to + tm + tp) / 3. También se podrá estimar un comportamiento uniforme para otros tipos de procesos, es decir que cualquier tiempo de duración tenga la misma probabilidad de ocurrencia en un rango determinado, y así se podrán identificar otros comportamientos.

Entonces lo que en verdad puede esperarse es que un cronograma de proyecto esté conformado por un conjunto de actividades entre las cuales unas actividades tengan un comportamiento PERT otras tengan un comportamiento lineal, algunas que tengan un comportamiento fijo o constante, etc.  Una buena opción para evaluar este cronograma es aplicando la simulación de Montecarlo que consiste en dar a cada actividad una duración dentro del rango de estimación y, con la secuencia que las interrelaciona, determinar una duración total del proyecto.

En estricto la simulación de Montecarlo aplicable a cualquier campo donde se requiere evaluar el riesgo tiene el siguiente proceso:

1. Diseñar el modelo lógico de decisión
2. Especificar distribuciones de probabilidad para las variables aleatorias relevantes.
3. Incluir posibles dependencias entre variables.
4. Muestrear valores de las variables aleatorias.
5. Calcular el resultado del modelo según los valores del muestreo (iteración) y registrar el resultado.
6. Repetir el proceso hasta tener una muestra estadísticamente representativa.
7. Obtener la distribución de frecuencias del resultado de las iteraciones.
8. Calcular media, desvío y curva de percentiles acumulados.

Las repeticiones de escenarios pueden ser de miles obteniendo tantas duraciones totales con las cuales se hace un cálculo estadístico que nos dirá la duración más frecuente, la que ocurre con una certeza de 100%, 90% 80% etc. y por lo tanto, si se tiene predefinido una duración T, podemos saber cuál es la certeza de que el cronograma cumpla con esta duración.

Las decisiones respecto a los resultados de la evaluación podrían ser las siguientes:
  • Si la duración total T calculada originalmente tiene una certeza de que se cumpla de 50%, hay que recalcular el cronograma o cambiar T por una duración con mayor certeza, por ejemplo mayor a 70%
  • Si la duración total T calculada originalmente tiene una certeza de que se cumpla de 90%, quizás puede uno sentirse más seguro.
  • etc.
Implementación en Excel

Debemos considerar que con cada cambio de las duraciones de las actividades interrelacionadas del proyecto se producirá una duración total. Resolver mediante fórmulas en celdas estos cálculos no es práctico más aún si debemos almacenar miles de duraciones para un postproceso, aquí es mejor recurrir a la programación mediante VBA y usar el excel como interfase para la introducción de datos y la presentación de resultados.

Un tipo de presentación puede ser el siguiente:


En una hoja deben indicarse las actividades y sus predecesoras y, una vez aplicado el proceso de calcular las duraciones para la cantidad de escenarios que se deseen, se puede visualizar los resultados con la siguiente presentación:


Teniendo definidas las formas en que se interrelacionan las actividades (secuenciado) se puede programar el cálculo de los tiempos de inicio y fin más tempranos y más tardíos de cada actividad y con esto encontrar la duración total del proyecto.

En el caso de cada actividad debe calcularse una duración aleatoria que cumpla con la característica que la defina. Es decir definir la distribución aleatoria que mejor represente su duración.

La distribución beta con los adecuados parámetros se ajusta a las probabilidades de duración, para el caso del método PERT, la distribución beta se presenta según lo siguiente:

Por ejemplo: sea una actividad con duración más optimista= 6, duración más pesimista= 16,

si la duración más probable es 13 la distribución beta adoptará la forma asimétrica a la derecha
si la duración más probable es 9 la distribución beta adoptará la forma asimétrica a la izquierda


Una vez determinados los parámetros para la distribución beta lo siguiente es ajustar para cada valor aleatorio entre 0 y 1 la duración de la actividad. Una propuesta de parámetros que produce la suficiente aproximación es:

p = 1 + 4.mo     q = 1 + 4.(1 - mo)

con    mo = (tm - to) / (tp - to)

En el Paper: Teaching Project Simulation in Excel Using PERT-Beta Distributions de Ron Davis, se expone unas ecuaciones que presentan una mejor aproximación:

        mo = (2 / 3 / (tp - to)) * (1 + 4 * (tm - to) * (tp - tm) / (tp - to) ^ 2)
        p = (tp + 4 * tm - 5 * to) * mo
        q = (5 * tp - 4 * tm - to) * mo

para el caso de la distribución triangular el ajuste de un numero aleatorio se puede hacer con las ecuaciones:

para a < (tm - to) / (tp - to)   d = to + raiz(a (tp - to) (tm - to))
para a >= (tm - to) / (tp - to)   d = tp - raiz((1- a) (tp - to) (tp - tm))

el caso en que una actividad tenga como duración cualquier tiempo entre dos valores mínimo y máximo es más simple y desde el número aleatorio a solo hay que aplicar:

d = (tmax - tmin).a + tmin  que corresponde a una distribución uniforme.

Un aplicativo que recoge lo explicado líneas arriba se puede descargar con el siguiente enlace:


algunas vistas de su presentación:



En la columna de distribuciones se hace click derecho para escoger el tipo de distribución:




Un cuadro de diálogo para cada tipo de distribución ayuda a introducir los parámetros:



Es todo por hoy.



Triangulación de una nube de puntos (una forma de calcularla) Parte 3

Posición de un punto respecto a un triángulo

Volviendo al asunto de determinar si un punto está dentro de un triangulo, ya sabemos que el área de los 3 triangulos que forman el punto con los lados del triangulo debe ser todos del mismo signo. El sentido antihorario de los vértices asegura el valor positivo.

Siendo estrictos diriamos que la suma de áreas de los 3 triángulos debe ser igual al área del triángulo que contiene al punto. Recalco que no es necesario calcular exactamente las áreas sino solo su signo.

En la figura siguiente se muestra las regiones en las que puede estar un punto respecto a un triángulo y nos sirve para entender lo que vamos a explicar:


Sea p el punto del cual queremos saber si está dentro del triangulo 0-1-2, la prueba requiere suponer que dicho punto está dentro del triangulo y por lo tanto se formarán los 3 triangulos siguientes:
  • p-0-1
  • p-1-2
  • p-2-0
Si el punto está en la región A, el triángulo p-2-0 será negativo y los otros 2 son positivos
Si el punto está en la región B, el triángulo p-0-1 será negativo
Si el punto está en la región C, el triángulo p-1-2 será negativo

Para comprobar lo dicho basta con verificar el cambio de sentido de los triangulos en cada una de las regiones, así se podrá ver que en las regiones intermedias AB, BC y CA son 2 los triángulos que cambian su sentido.

Por lo tanto, para los efectos que nos demanda el algoritmo es posible hacer el código de manera que podemos no solo saber si un punto está dentro de un triangulo sino que podemos saber de que lado del triangulo se encuentra.

Entonces esto justifica sostener que el registro de cada triangulo debe contener los vértices que lo conforman y tambien los triangulos que son vecinos a él. Por ejemplo en la siguiente figura se muestra con fondo ámbar la numeración de triángulos y con fondo transpartente la numeración de vértices.



Los indices que numeran los vértices son los que corresponden al arreglo o lista de puntos con sus coordenadas. Para los triangulos adoptamos la siguiente forma de describirlo: (para el triangulo 6)
           ( (10 9 8) 8 9 7)

esta es una lista de 4 elementos, el primer elemento: (10 9 8) es el triángulo principal descrito con los vértices 10, 9 y 8, y los otros tres números representan los indices de los triángulos vecinos ordenados de modo que, para el ejemplo: el triangulo 8 es opuesto al vértice 10, el triangulo 9 es opuesto al vértice 9 y el triángulo 7 es opuesto al vértice 8. Estos tres últimos triángulos sons vecinos y obviamente tienen como arista común los lados 8-9, 10-8 y 9-10 respectivamente.

De acuerdo a la forma adoptada para referirnos a cada triangulo, resulta mas sencillo extraer datos del triangulo principal desde CAR y los identificadores del "vecindario" desde CDR. El tipo de registro presentado es de 4 triangulos.

Con este escenario es posible escribir el código para mejorar el algoritmo de triangulación en los puntos más cruciales:

- hallar el triangulo que contiene un nuevo punto lo más rapido posible

- efectuar la legalización de los nuevos triángulos determinando inmediatamente los triángulos "vecinos" o los que tienen una arista común.

La siguiente figura ilustra el proceso que sigue la busqueda del triángulo que contiene a un punto p que se suministra al código.



Si el proceso de verificación se inicia en el triangulo 0, la verificación indica que el siguiente triángulo a revisar es el 1, la verificación en este indicará que debe evaluarse el triangulo 2, y así hasta llegar al triángulo de color verde en donde al verificar se encontrará que es el triangulo buscado.

Una vez hallado el triangulo contenedor se descompone en 3 triángulos "hijos" y se actualiza su relación con el "vecindario" del triangulo original. Esto quiere decir que debe registrarse en los triángulos vecinos a los nuevos hijos ("herederos" del espacio del triangulo original). Algo similar ocurre cuando se hace el "flip" que legaliza 2 triangulos, se produce un cambio en la relación de los "vecinos" con los triangulos "flipeados" que debe quedar registrado.

El mecanismo detallado para automatizar estas operaciones es algo engorroso de explicar. Es muy importante mantener algún patrón de numeración de vértices, hecho que facilitará algunas operaciones.

Utilidad

La triangulación de una nube de puntos tiene una aplicación práctica en la ingeniería civil para la determinación de un modelo digital del terreno (MDT). Los puntos del terreno resultan del levantamiento topográfico o de algun otro proceso. El MDT es la base sobre la que se automatizan algunas operaciones que asisten al diseño de ingeniería.

Hace no mucho compartí en algunos foros especializados un compilado de utilidades para operaciones de topografía en gabinete que contiene como elemento principal una implementación de la triangulación de la nube de puntos tal como fue expuesta aquí, una versión recompilada puede obtenerse en el siguiente enlace:



una vez cargado el archivo en el Autocad se tendrán disponibles unos comandos entre los que se destaca los siguientes:
  • importar registros de puntos de txt
  • generar la triangulación como 3dface y lineas tin
  • dibujo de las curvas de nivel
algunas de sus interfaces se presentan con las siguientes vistas:

Este cuadro se muestra ejecutando el comando PUNTOS


Este cuadro se muestra ejecutando el comando REUBIK

Este cuadro se muestra ejecutando el comando CONF-MDT


Una vez cargados los puntos se utiliza el comando TIN que pedirá seleccionar en pantalla el conjunto de puntos a triangular.

las curvas de nivel se dibujan como lineas rectas para cada triángulo, cuando sea necesario se puede usar el comando SUAVIZA para obtener un suavizado de las curvas de nivel, sobre las curvas maestras se escribirán las cotas en que se encuentran.

Con la triangulación ya es posible aplicar herramientas propias del diseño en ingeniería. Donde más se aprovecha es en el diseño vial de carreteras o caminos ya que se facilitan los cálculos de movimientos de tierras entre otras operaciones.

Tambien trataremos el desarrollado de un conjunto de herramientas para el diseño geométrico.


Triangulación de una nube de puntos (una forma de calcularla) Parte 2

Aplicando Autolisp...

La siguiente funcion calcula el área de un triángulo proporcionándole sus tres vértices, (en estricto el área es la mitad del valor que devuelve esta función):

(defun signarea (pa pb pc)
   (+ (* (car pa) (- (cadr pb) (cadr pc)))
      (* (car pb) (- (cadr pc) (cadr pa)))
      (* (car pc) (- (cadr pa) (cadr pb))))
)

Una Implementación inicial del algoritmo:

Un algoritmo de triangulación consistente en encontrar el triángulo que contiene cada nuevo punto y, una vez hallado, formar tres nuevos triángulos, los que deben pasarse por la verificación de la condición de Delaunay, se presenta en la siguiente forma:

la variable: listadepuntos es la que contiene a la nube de puntos.

(defun tin ()
  ;;(obtener_xmax_xmin_ymax_ymin)
ß para no distraernos de lo principal…
  ;; lo siguiente determina el triángulo mayor
  ;; o sea el que contiene a todos los puntos
  (setq pa (list xmin (+ ymax (* 2.0 (- ymax ymin))) 0)
        pb (list (- xmin (* 2.5 (- xmax xmin)))
                 (- ymin (* 2.5 (- ymax ymin))) 0)
        pc (list (+ xmax (* 2.5 (- xmax xmin))) ymin 0) )
  (setq listadepuntos (append (list pa pb pc) listadepuntos))

;; primeros 3 triángulos
  (setq tr1 '(3 0 1) tr2 '(3 1 2) tr3 '(3 2 0) )
  (setq grupo_triang (list tr1 tr2 tr3))
  (setq ii 4
        uu (length listadepuntos))
  (while (< ii uu)
    (setq eli (halla_tri))
    (setq grupo_triang (eliminar eli grupo_triang))
    (setq tr1 (list ii (car eli) (cadr eli))
          tr2 (list ii (cadr eli) (caddr eli))
          tr3 (list ii (caddr eli) (car eli)) )
    (setq gruponuevo (list tr1 tr2 tr3))
    (legalizar gruponuevo)
;; la siguiente línea añade los nuevos triángulos legalizados a
;; la triangulación. En el proceso de legalización se forman
;; nuevos triángulos en lugar de otros por la acción del flip.
    (setq grupo_triang (append grupo_triang gruponuevo))
    (eliminar el triangulo que contiene el punto nuevo)
    (setq ii (+ ii 1))
    )
  (dibuja_red_tin)
)

;; lo siguiente devuelve el triángulo que contiene al nuevo punto
(defun halla_tri ( / tr00 jj encontrado)
  (setq jj 0

  (while (not encontrado)
    (setq tr00 (nth jj grupo_triang)
          encontrado
           (and (> (signarea (list ii (car tr00) (cadr tr00)))  0.0)
                (> (signarea (list ii (cadr tr00) (caddr tr00))) 0.0)
                (> (signarea (list ii (caddr tr00) (car tr00)))) 0.0) )
          jj (1+ jj)

       )
  )
  tr00


(defun signarea (a b c)

  (setq pa (nth a listadepuntos)
        pb (nth b listadepuntos)
        pc (nth c listadepuntos)
   )
  (+ (* (car pa) (- (cadr pb) (cadr pc)))
     (* (car pb) (- (cadr pc) (cadr pa)))
     (* (car pc) (- (cadr pa) (cadr pb))) ) 
 )

El código mostrado  lineas arriba complementado con las operaciones de legalización, eliminación de cada triangulo que contiene un punto nuevo y dibujado, nos produce la triangulación buscada.

El procedimiento que se aplica corresponde al algoritmo incremental aleatorio, en primer lugar se define un triangulo grande que contiene a toda la nube y luego se van añadiendo los puntos uno por uno, en cada incorporación de puntos se verifica que los triangulos formados cumplan con la condición Delaunay. Al terminar se hace el dibujado de los lados de la triangulación (líneas TIN) cuidando que no sean las que conectan los vértices del supertriángulo y que no se repitan.

De todas las operaciones a efectuar con este código; la de encontrar el triángulo que contiene a un nuevo punto se convierte en crucial si se quieren procesar miles de puntos. La función de búsqueda presentada líneas arriba es demasiado ordinaria. Es necesario implementar funciones más eficientes, una mejor opción es aprovechar una de las funciones de Visual Lisp como lo siguiente:

(defun halla_tri ( / x ) 
  (car (vl-member-if '(lambda (x)
              (and (>= (signarea (list ii (car x) (cadr x)))  0.0)
                   (>= (signarea (list ii (cadr x) (caddr x))) 0.0)
                   (>= (signarea (list ii (caddr x) (car x)))) 0.0) ) )
           
           grupo_triang) )

 )

A esta solución podemos añadir un preordenamiento de la lista de puntos según las coordenadas X, tambien ayuda que los triangulos nuevos se añadan al comienzo de la lista: 

(setq grupo_triang (append gruponuevo grupo_triang) )

Sobre el "legalizado"

Una triangulación será correcta si se cumple con la condición de Delaunay, entonces cómo lo verificamos?.
Si los vértices de cada triángulo conforman una circunferencia dentro de la cual no se encuentra ningún otro punto de la nube, se dirá que se trata de una triamgulación Delaunay.
Conozco 2 maneras de efectuar la verificación:

1. Hallar el centro de la circunferencia circunscrita, por consiguiente el radio de dicha circunferencia y verificar la distancia del punto a evaluar hasta el circuncentro, si la distancia es menor que el radio, el punto está dentro de la circunferencia y debe efectuarse la modificación de la triangulación.

2. Resolver la determinante haciendo la verificación:

\begin{vmatrix}
A_x & A_y & A_x^2 + A_y^2 & 1\\
B_x & B_y & B_x^2 + B_y^2 & 1\\
C_x & C_y & C_x^2 + C_y^2 & 1\\
D_x & D_y & D_x^2 + D_y^2 & 1
\end{vmatrix} = \begin{vmatrix}
A_x - D_x & A_y - D_y & (A_x - D_x)^2 + (A_y - D_y)^2 \\
B_x - D_x & B_y - D_y & (B_x - D_x)^2 + (B_y - D_y)^2 \\
C_x - D_x & C_y - D_y & (C_x - D_x)^2 + (C_y - D_y)^2 \\
\end{vmatrix} > 0

donde A, B, C son los vértices del triángulo y D es el punto que se desea verificar. Si resulta mayor que cero, el punto D está dentro de la circunferencia

También aquí el procedimiento debe ser eficiente. Si consideramos que no conocemos de antemano los puntos candidatos a verificar por cada triangulo, tendremos que verificar todos los lados de la triangulación lograda hasta encontrar la del triangulo que tiene el lado comun. Conforme crece la triangulación el proceso se hace más lento y se agrava si la cantidad de puntos de la nube es grande (> 2000).

Una versión del algoritmo en Autolisp se puede descargar de aqui:

triangulacion en lisp

Esta versión contiene un mejoramiento de los procedimientos y permite obtener triangulaciones más rápidas que con los procedimientos no refinados mostrados líneas arriba.

Es posible plantearse procedimientos más elaborados que mejoren la eficiencia del algoritmo. Se debe considerar la necesidad de que por cada triángulo se almacenen más registros que los 3 vértices que lo conforman. Mas adelante explicaré un planteamiento alternativo:

Cada registro de triángulo tiene los vértices que lo conforman y los "vecinos" por cada lado.

(continuará...)