$TITLE 'Optimizacion bajo incertidumbre aplicado a la industria del petroleo'
$OFFUPPER
$ONEMPTY
$OFFLISTING
$OFFSYMXREF
$OFFSYMLIST
OPTION SOLPRINT = ON;
SETS
 P   'Productos' / Cl 'Crudo ligero', Cm 'Crudo medio',
   Cp 'Crudo pesado', Gl  'Gasolina', Ke  'Queroseno', Go  'Gasoleo' /
 Po   'Productos que pueden transportar los oleoductos' / Cl, Cm, Cp /
 I   'Centros: extraccion crudo, almacenamiento, refino, puertos' /
   E1 'Pozos de petroleo area E1', E2 'Pozos de petroleo area E2',
   E3 'Pozos de petroleo area E3',
   N1 'Estacion de bombeo e interconexion de oleoductos N1',
   R1 'Refineria R1 en estudio', B1 'Puerto y centro de almacenamiento B1',
   C1 'Puerto y centro de venta de productos y crudos'  /
 M   'Modos de funcionar una refineria' /
   Ml 'Refino de crudo ligero Cl', Mm 'Refino de crudo medio Cm'
   Mp 'Refino de crudo pesado Cp' /
 Ro 'Tramos de oleoducto' / N1B1, E1N1, E2N1, E2B1, E2R1, E2E3, E3R1, B1R1 /
 F   'Sentido de flujo de oleoductos'  / D 'Directo', I 'Inverso'  /
 RIo(ro,f,i) 'Centros de origen de flujo en tramos oleoducto' /
   N1B1 . D . N1, E1N1 . D . E1, E2N1 . D . E2, E2B1 . D . E2,
   E2R1 . D . E2, E2E3 . D . E2, E3R1 . D . E3, B1R1 . D . B1
   N1B1 . I . B1, E1N1 . I . N1, E2N1 . I . N1, E2B1 . I . B1,
   E2R1 . I . R1, E2E3 . I . E3, E3R1 . I . R1, B1R1 . I . R1  /
 RFo(ro,f,i) 'Centros de destino de flujo en tramos oleoducto' /
   N1B1 . I . N1, E1N1 . I . E1, E2N1 . I . E2, E2B1 . I . E2,
   E2R1 . I . E2, E2E3 . I . E2, E3R1 . I . E3, B1R1 . I . B1
   N1B1 . D . B1, E1N1 . D . N1, E2N1 . D . N1, E2B1 . D . B1,
   E2R1 . D . R1, E2E3 . D . E3, E3R1 . D . R1, B1R1 . D . R1  /
 B   'Buque tanque de productos y crudos' /
   Bc 'Buques de crudo de 250.000 TPM', Bp 'Buques de productos de 70.000 TPM'/
 Rb  'Trayectos de buques'  /
   B1C1 'Trayecto de buque de B1 a C1', R1C1 'Trayecto de buque de R1 a C1' /
 RIb(rb,i) 'Centros origen flujo en trayectos buque' / B1C1 . B1, R1C1 . R1 /
 RFb(rb,i) 'Centros destino flujo en trayectos buque' / B1C1 . C1, R1C1 . C1 /
 Pb(b, p)  'Productos que puede transportar cada buque' /
    Bc . (Cl, Cm, Cp), Bp . (Gl, Ke, Go)   /
 G   'Escenarios'  / G1 'Escenario 1', G2 'Escenario 2' /    ;
ALIAS (i, j), (p, p1), (b, b1) ;
PARAMETER
 cFe(i) 'Coste fijo nuevos pozos extraccion crudo' / E3 70  /
 qE_X(g,i,p) 'Nueva capacidad extraccion crudo' /
 G1 . E3 . Cl   30 , G2 . E3 . Cl   40  / ;
TABLE
 qE_A(g,i,p) 'Capacidad actual extraccion crudo'
           Cm     Cp
 G1 . E1          10.
 G2 . E1          6.
 G1 . E2   20.
 G2 . E2   35.        ;
TABLE
 cVe(g,i,p)  'Coste variable de extraccion de crudo'
           Cl     Cm     Cp
 G1 . E1                 43.5
 G2 . E1                 65.8
 G1 . E2          43.1
 G2 . E2          54.2
 G1 . E3   54.1
 G2 . E3   54.8                ;
*  Refino
TABLE
etaR(i,m,p) 'Producto p (>0, producido; <0, consumido) en refineria R'
         Cl     Cm     Cp    Gl     Ke     Go
R1 . Ml  -1.                 0.59   0.14   0.25
R1 . Mm         -1.          0.44   0.19   0.35
R1 . Mp                -1.   0.25   0.18   0.55  ;
TABLE
cVr(i,m)   'Coste variable de refino'
    Ml     Mm     Mp
R1  25     27     31 ;
PARAMETER
cFr(i)  'Coste fijo de la refineria mas el terminal' /
  R1   84.  /
qR_X(i) 'Capacidad maxima de refino' / R1   20. / ;
TABLE
 cVo(ro,p)   'Coste variable de transporte por oleoducto'
          Cl     Cm     Cp
    N1B1  3.32   3.34   3.37
    E1N1  4.45   4.46   4.48
    E2N1  2.25   2.26   2.27
    E2B1  4.44   4.45   4.47
    E2R1  3.3    3.31   3.33
    E2E3  3.36   3.37   3.38
    E3R1  2.2    2.21   2.22
    B1R1  6.28   5.79   6.29    ;
PARAMETER
 qO_A(ro)  'Caudal maximo actual oleoductos' /
    N1B1  40, E1N1  20,  E2N1  30  /
 qO_X(ro)  'Caudal maximo nuevos oleoductos' /
    E2B1  30, E2R1  50, E2E3  50, E3R1  60, B1R1  30 /  ;
PARAMETER
 cFo(ro)     'Coste fijo nuevos oleoductos'   /
    E2B1  61, E2R1  45, E2E3  55, E3R1  34, B1R1  55  / ;
*  Buques
PARAMETER
cFb(g, b)      'Coste fijo buque' /
  G1 . Bc 24.5, G1 . Bp 24.7, G2 . Bc 22.3, G2 . Bp 22.6  /
cVb(g, b, rb)  'Coste variable viaje redondo buque' /
  G1 . Bc . B1C1  4.45, G1 . Bc . R1C1  4.42, G1 . Bp . B1C1  3.23,
  G1 . Bp . R1C1  3.21, G2 . Bc . B1C1  5.51, G2 . Bc . R1C1  4.47,
  G2 . Bp . B1C1  3.30, G2 . Bp . R1C1  3.27    /
hVr(b, rb)  'Duracion viaje redondo buque' /
  Bc . B1C1  210, Bc . R1C1  208, Bp . B1C1  205, Bp . R1C1  201  /
qB(b)      'Carga util maxima del buque' /
  Bc 0.240, Bp 0.065   /
tuaB(b)      'Tiempo neto util anual de un buque' /
  Bc 8300,  Bp 7900  /  ;
*  Comercializacion
TABLE
 pV(g,i,p)   'Precio de venta de p en i para escenario g'
           Cl     Cm     Cp     Gl     Ke     Go
 G1 . R1   075.   060.   045.   120.   110.   140.
 G2 . R1   125.   130.   125.   160.   175.   170.
 G1 . B1   080.   075.   060.
 G2 . B1   140.   110.   190.
 G1 . C1   090.   105.   088.   155.   181.   138.
 G2 . C1   165.   115.   116.   194.   193.   188. ;
TABLE
 qD_X(g,i,p) 'Demanda maxima p en i para escenario g'
           Cl     Cm     Cp     Gl     Ke     Go
 G1 . R1   5.     3.     4.     1.5    2.1    4.8
 G2 . R1   7.     8.     6.     2.     3.     6.
 G1 . B1   10.    13.    11.
 G2 . B1   12.    16.    14.
 G1 . C1   8.     9.     7.     5.3    6.     6.
 G2 . C1   7.     6.     8.6    3.5    8      7.     ;
PARAMETER
 wG(g) 'Probabilidad de que suceda el escenario g' / G1 0.75 , G2  0.25 / ;
SET
PoNo(p) 'Productos que no pueden circular por oleoductos'
Gid(g,i,p) 'Centros con demanda de p en escenario g'
Gie(g,i,p) 'Centro con capacidad de extraccion de crudo'
IeN(i) 'Centro con capacidad nueva de extraccion de crudo'
Gir(i,m) 'Centro con capacidad de refino modo m'
Gi(i) 'Centro con capacidad de refino modo m'  ;
PoNo(p) = NOT Po(P) ;
IeN(i) = SUM((g,p)$(qE_X(g,i,p) gt EPS), YES) ;
Gid(g,i,p) = YES$(qD_X(g,i,p) gt EPS) ;
Gie(g,i,p) = YES$((qE_X(g,i,p)+qE_A(g,i,p)) gt EPS) ;
Gir(i,m) = YES$(qR_X(i)$(SUM(p,ABS(etaR(i,m,p))) gt EPS)) ;
Gi(i) = SUM(m$Gir(i,m), YES) ;
* El modelo
FREE VARIABLES
Benef    'Beneficios anuales (millones E)'  ;
POSITIVE VARIABLES
Fo_grfp(g,ro,f,p) 'Flujo en oleoducto r, f, p para escenario g (millones t/a)'
Fb_grp(g,b,rb,p)  'Flujo por buque r, f, p para escenario g (millones t/a)'
E_gip(g,i,p)     'Crudo p extraido en centro i en escenario g (millones t/a)'
R_gim(g,i,m)     'Crudo refinado en centro i en escenario g (millones t/a)'
D_gip(g,i,p)     'Demanda de (g, i, p) satisfecha (millones t/a)'
Nv(g,b,rb)       'Numero de viajes de buque b en trayecto r'
Nb(g,b)          'Numero de buques necesarios de tipo b'
BINARY VARIABLES
Yo(ro)    '(0=no, 1=si) construimos el oleoducto r'
Ye(i)     '(0=no, 1=si) construimos los pozos de extraccion del area i'
Yr(i)     '(0=no, 1=si) construimos la refineria del area i'
EQUATIONS
Fobj        'Funcion objetivo'
Ec1(g,i,p)  'Balance producto p en centro i para g'
Ec2(g,i,p)  'Capacidad de extraccion de crudo'
Ec3(g,i)    'Capacidad de refino'
Ec4(g,ro)   'Capacidad oleoductos'
Ec5(g,b,rb) 'Viajes de buques por trayecto'
Ec6(g,b)    'Numero de buques necesarios'    ;
Fobj.. Benef - SUM((g,i, p)$Gid(g,i,p), wG(g)*pV(g,i,p)*D_gip(g,i,p))
       + SUM(i$IeN(i), cFe(i)*Ye(i))
       + SUM(i$(qR_X(i) gt EPS), cFr(i)*Yr(i))
       + SUM((g,i,p)$Gie(g,i,p), wG(g)*cVe(g,i,p)*E_gip(g,i,p))
       + SUM((g,i,m)$Gir(i,m), wG(g)*cVr(i,m)*R_gim(g,i,m))
       + SUM(ro$(qO_X(ro) gt EPS), cFo(ro)*Yo(ro))
       + SUM((g,ro,f,p), wG(g)*cVo(ro,p)*Fo_grfp(g,ro,f,p))
       + SUM((g,b,rb), wG(g)*cVb(g, b, rb)*Nv(g,b,rb))
       + SUM((b, g), wG(g)*cFb(g, b)* Nb(g,b)) =E= 0 ;
Ec1(g,i,p).. E_gip(g,i,p)$Gie(g,i,p) - D_gip(g,i,p)$Gid(g,i,p)
       + SUM(m$(ABS(etaR(i,m,p)) gt EPS), etaR(i,m,p)*R_gim(g,i,m))
       - SUM((ro,f)$RIo(ro,f,i), Fo_grfp(g,ro,f,p))
       + SUM((ro,f)$RFo(ro,f,i), Fo_grfp(g,ro,f,p))
       - SUM((rb,b)$(RIb(rb,i)$Pb(b,p)), Fb_grp(g,b,rb,p))
       + SUM((rb,b)$(RFb(rb,i)$Pb(b,p)), Fb_grp(g,b,rb,p)) =E= 0 ;
Ec2(g,i,p)$Gie(g,i,p).. E_gip(g,i,p) - qE_X(g,i,p)*Ye(i) =L= qE_A(g,i,p) ;
Ec3(g,i)$Gi(i).. SUM(m$Gir(i,m), R_gim(g,i,m)) - qR_X(i)*Yr(i) =L= 0 ;
Ec4(g,ro).. SUM((p,f), Fo_grfp(g,ro,f,p))
        - (qO_X(ro)*Yo(ro))$(qO_X(ro) gt EPS) =L= qO_A(ro) ;
Ec5(g,b,rb).. SUM(p$Pb(b,p), Fb_grp(g,b,rb,p)) - qB(b)* Nv(g,b,rb) =L= 0;
Ec6(g,b).. SUM(rb, hVr(b, rb)*Nv(g,b,rb)) - tuaB(b)*Nb(g,b) =L= 0;
* Las opciones
  OPTION ITERLIM = 4000000 ;
  OPTION LIMCOL  = 0 ;
  OPTION LIMROW  = 0 ;
  OPTION OPTCR   = 0.0000001 ;
  OPTION RESLIM  = 360000 ;
  OPTION SOLPRINT = ON ;
* elegir una de las dos opciones
  OPTION MIP   = CPLEX  ;
*  OPTION MIP   = XPRESS  ;
* Elaboracion del modelo
MODEL Petroleo /ALL/ ;
* Limites cantidades vendidas
E_gip.UP(g,i,p) = qE_A(g,i,p) + qE_X(g,i,p) ;
D_gip.UP(g,i,p) = qD_X(g,i,p) ;
Fo_grfp.FX(g,ro,f,p)$PoNo(p) = 0 ;
* optimizacion
SOLVE Petroleo MAXIMIZING Benef USING MIP ;
* salida de resultados
DISPLAY Benef.L, Fo_grfp.L, Fb_grp.L, E_gip.L, R_gim.L ;
DISPLAY D_gip.L, Nv.L, Nb.L, Yo.L, Ye.L, Yr.L ;
FILE Arc_SalRes 'Archivo salida de resultados' /Resultados.txt/;
PUT Arc_SalRes ;
Arc_SalRes.pc = 5 ;
 PUT 'Beneficios anuales', Benef.L :16:5 / ;
 PUT 'Nuevo centro extraccion crudo' / ;
 LOOP(i$(Ye.L(i) gt EPS),
   PUT i.tl, Ye.L(i) / ;
 ) ;
 PUT 'Nuevas refinerias' / ;
 LOOP(i$(Yr.L(i) gt EPS),
   PUT i.tl, Yr.L(i) / ;
 ) ;
 PUT 'Nuevo tramo oleoducto' / ;
 LOOP(ro$(Yo.L(ro) gt EPS),
   PUT ro.tl, Yo.L(ro) / ;
 ) ;
PARAMETER aux2(ro);
aux2(ro) = qO_A(ro)+qO_X(ro)*Yo.L(ro) ;
 PUT 'Escenario','Tramo','Sentido','Producto','Flujo oleoducto','Max flujo'/ ;
 LOOP((g,ro,f,p)$(Fo_grfp.L(g,ro,f,p) gt EPS),
   PUT g.tl, ro.tl, f.tl, p.tl, Fo_grfp.L(g,ro,f,p) :16:5, aux2(ro) :16:5 / ;
 ) ;
 PUT 'Escenario', 'Buque', 'Trayecto', 'Producto', 'Flujo buque' / ;
 LOOP((g,b,rb,p)$(Fb_grp.L(g,b,rb,p) gt EPS),
   PUT g.tl, b.tl, rb.tl, p.tl, Fb_grp.L(g,b,rb,p) :16:5 / ;
 ) ;
PARAMETER Aux(g,i,p);
Aux(g,i,p) = qE_X(g,i,p)*Ye.L(i)+qE_A(g,i,p);
 PUT 'Escenario', 'Pozo petr', 'Crudo', 'Flujo extraido', 'Max capacidad' / ;
 LOOP((g,i,p)$(E_gip.L(g,i,p) gt EPS),
   PUT g.tl, i.tl, p.tl, E_gip.L(g,i,p) :16:5, Aux(g,i,p) :16:5 / ;
 ) ;
 PUT 'Escenario', 'Refineria', 'Modo ref', 'Nivel refino', 'Max capacidad' / ;
 LOOP((g,i,m)$(R_gim.L(g,i,m) gt EPS),
   PUT g.tl, i.tl, m.tl, R_gim.L(g,i,m) :16:5, qR_X(i) :16:5 / ;
 ) ;
 PUT 'Escenario', 'Centro', 'Producto', 'Demanda serv', 'Max demanda' / ;
 LOOP((g,i,p)$(D_gip.L(g,i,p) gt EPS),
   PUT g.tl, i.tl, p.tl, D_gip.L(g,i,p) :16:5, qD_X(g,i,p) :16:5 / ;
 ) ;
 PUT 'Escenario', 'Buque', 'Trayecto', 'Num viajes redondos' / ;
 LOOP((g,b,rb)$(Nv.L(g,b,rb) gt EPS),
   PUT g.tl, b.tl, rb.tl, Nv.L(g,b,rb) :16:5 / ;
 ) ;
 PUT 'Escenario', 'Buque', 'Num buques' / ;
 LOOP((g,b)$(Nb.L(g,b) gt EPS),
   PUT g.tl, b.tl, Nb.L(g,b) :16:5 / ;
 ) ;
PUTCLOSE Arc_SalRes ;
* fin del programa