korpet.gms : Investment Planning in the Korean Oil-Petro Industry

Description

This problem addresses two questions in the Korean oil refining and
petrochemical industries: a) the optimal time for capacity expansion,
and b) options available to deal with increasingly strict government
anti-pollution regulations. An interesting problem is that due to
inflexibility in production scheduling (as stated) in the petrochemical
sector it is not possible to satisfy market requirements without
producing excess amounts of certain products which then have to be
disposed of.


Large Model of Type : MIP


Category : GAMS Model library


Main file : korpet.gms

$title Korea: Dynamic Oil-Petro Model (KORPET,SEQ=48)

$onText
This problem addresses two questions in the Korean oil refining and
petrochemical industries: a) the optimal time for capacity expansion,
and b) options available to deal with increasingly strict government
anti-pollution regulations. An interesting problem is that due to
inflexibility in production scheduling (as stated) in the petrochemical
sector it is not possible to satisfy market requirements without
producing excess amounts of certain products which then have to be
disposed of.


Suh, J S, An Investment Planning Model for the Oil-Refining and
Petro-chemical Industries in Korea. Tech. rep., Center for Economic
Research, University of Texas, 1981.

Keywords: mixed integer linear programming, oil refining, investment planning,
          energy economics
$offText

$sTitle Set Definitions
Set
   i 'plant locations'
     / ulsan   'kyungsangnamdo'
       yosu    'chollanamdo'
       inchon  'kyunggido'      /
   j 'demand regions'
     / seoul   'capital city'
       inchon  'kyunggido'
       kwangju 'chollanamdo'
       taegu   'kyungsangbukdo'
       pusan   'kyungsangnamdo'
       ulsan   'kyungsangnamdo'
       yosu    'chollanamdo'    /
   m 'productive units'
     / ad      'atmospheric distillation unit'
       cr      'catalytic reformer'
       cc      'catalytic cracker'
       ds      'desulfurizer'
       sc      'steam cracker'
       bx      'butadiene extractor'
       pf      'platfiner'
       ar      'aromatics unit'
       hd      'hydrodealkylator'              /
   p 'processes'
     / crad    'crude oil atmospheric distillating process'
       nacr    'naphtha reforming process'
       focc    'fuel oil catalytic cracking process'
       rsds    'straight run residuum desulfurizing process'
       nasc    'naphtha steam cracking process'
       btbx    'butadine extraction process'
       pgpf    'raw pyrolysis gas platfining process'
       amna    'aromatics processing using reformed naphtha'
       amtp    'aromatics processing using treated pyrolysis gas'
       tohd    'toluene hydrodealkylating process'                /
   c 'commodities'
     / sa      'saudi arabian crude'
       adga    'gasoline from atmospheric distillation'
       adna    'naphtha from atmospheric distillation'
       adfo    'fuel oil from atmospheric distillation'
       adke    'kerosene from atmospheric distillation'
       adrs    'residuum from atmospheric distillation'
       crna    'catalytic reformed naphtha'
       crga    'catalytic reformed gasoline'
       ccga    'catalytic cracked gasoline'
       ccfo    'catalytic cracked fuel oil'
       dsrs    'desulfurized residuum'
       scbt    'steam cracked raw c4'
       scrp    'raw pyrolysis gas from steam cracking'
       pftp    'treated pyrolysis gas  from platfining'
       pg      'premium gasoline'
       rg      'regular gasoline'
       di      'distillate'
       rf      'residual fuel oil'
       et      'ethylene'
       pp      'propylene'
       bt      'butadiene'
       be      'benzene'
       to      'toluene'
       mx      'mixed xylene'                             /
   cfr(c)    'final products from refineries'             / pg, rg, di, rf /
   cfrerf(c) 'refinery products except residual fuel oil' / pg, rg, di     /
   cfrrf(c)  'residual fuel oil only'                     / rf /
   cfp(c)    'final products from petrochemical plants'   / et, pp, bt, be, to, mx /
   ci(c)     'intermediate products'
             / adga, adna, adfo, adke, adrs, crna, crga
               ccga, ccfo, dsrs, scbt, scrp, pftp       /
   cs(c)     'interplants shipments' / adna, crna /
   cr(c)     'raw materials'         / sa         /
   cb(c)     'commodities for blending refinery products'
             / adga, adfo, adke, adrs, crga
               crna, ccga, ccfo, dsrs        /
   cfcb(c,c) 'allowed blending combinations'
             /(pg,rg).(adga,adna,crga,crna,ccga)
               di    .(adke,adfo,ccfo)
               rf    .(adfo,ccfo,adrs,dsrs)     /
   crcr(c,c) 'crude oil combination in processing' / sa.sa /
   q         'quality specification attributes'
             / oc  'research octane number'
               va  'vapor pressure'
               su  'sulfur content (weight percentage)' /
   t         'time periods'
             / 1980-84, 1985-89, 1990-94, 1995-99 /
   te(t)     'expansion  period'
             / 1985-89, 1990-94, 1995-99 /
   cf(c)     'final products from refineries and petrochemical plants'
   crcf(c)   'raw materials and final products';

Alias (t,tp), (i,ip), (te,tep);

cf(cfr)  = yes;
cf(cfp)  = yes;
crcf(cr) = yes;
crcf(cf) = yes;
display cf, crcf;

$sTitle Data
Scalar year   'years per time period' / 5 /;

Parameter
   midyear(t) 'year which is in the middle of period'
   ts(t,t)    'time summation matrix';

midyear(t) = 1977 + year*ord(t);
ts(t,tp)$(ord(tp) <= ord(t)) = 1;
display ts;

Table a(c,c,p) 'input-output coefficients'
                crad   nacr   focc   rsds
   sa.sa       -1.00
   sa.adga       .04
   sa.adna       .14  -1.00
   sa.adke       .09
   sa.adfo       .18         -1.00
   sa.adrs       .54                -1.00
   sa.crga              .16
   sa.crna              .56
   sa.ccga                     .51
   sa.ccfo                     .25
   sa.dsrs                           1.00

   +            nasc   btbx   pgpf   amna   amtp   tohd
   sa.adna     -1.00
   sa.scbt       .08  -1.00
   sa.scrp       .23         -1.00
   sa.crna                          -1.00
   sa.pftp                     .68         -1.00
   sa.et         .29
   sa.pp         .15
   sa.bt                .53
   sa.to                              .26    .24  -1.00
   sa.be                              .40    .39    .63
   sa.mx                              .15    .16       ;

Table b(m,p) 'capacity utilization coefficients'
              crad nacr focc rsds nasc btbx pgpf amna amtp tohd
   ad          1.0
   cr               1.0
   cc                    1.0
   ds                         1.0
   sc                              1.0
   bx                                   1.0
   pf                                        1.0
   ar                                             1.0  1.0
   hd                                                       1.0;

Table ka(m,i) 'initial capacity (1000 tpy)'
               inchon  ulsan   yosu

   ad           3702.  12910.  9875.
   cr            369.   1200.   790.
   cc
   ds
   sc                    517.  1207.
   bx                     45.    94.
   pf                    201.
   ar                    181.   148.
   hd                    180.       ;

Parameter
   rket1(i) 'rate of capacity expansion at time period 1' / inchon .3, ulsan .4, yosu .6 /
   kat1     'planned capacity in time period 1';

kat1(m,i) = ka(m,i)*(1 + rket1(i));
display kat1;

Parameter op(p) 'operating cost      (us $ per ton)'
                / crad    1.55, nacr    9.71, focc  3.5 , rsds   14
                  btbx  101.82, nasc  143.04, pgpf 50.25, amna  190.87
                  amtp  190.87, tohd  645.99                           /
* operating cost includes catalyst cost,chemical cost,labor
* cost,depreciation cost,utility costs (steam,electricity,
* cooling water) and maintenance cost.

Table d(c,j) 'ratio of regional demand (percentage)'
* all petrochemical final products are demanded by the downstream petrochemical plants which are
* located very cloase to the refineries and petrochemical plants.  the demand pattern between regions
* for final products is assumed to remain constant over time.
              inchon  seoul
   pg            1.7   76.4
   rg            3.6   53.7
   di           22.4   31.4
   rf           68.0

   +         kwangju  taegu  pusan  ulsan  yosu
   pg           2.3    10.6    9.0
   rg           10.4   19.6   12.7
   di           8.1    19.4   18.7
   rf           5.0    15.8    9.6     .9    .7
   et                                38.1  61.9
   pp                                30.2  69.8
   bt                                28.6  71.4
   be                                37.8  62.2
   to                                30.8  69.2
   mx                                60.9  39.1;

Parameter
   tqcf(c)  'total quantity demanded     (1000 tpy)'
            / pg  335.9, rg  783.1, di  5391.2, rf  12696.1
              et  184.5, pp  219.3, bt    35.6, be    110.5
              to   84.1, mx   59.8                          /
* the growth rate of demand for pg,rg and petrochemical products
* is assumed to be 15 per cent per year and the growth rate of
* demand for the rest of the refinery products is assumed to be
* 10 per cent per year.

   grfp(c)  'growth rate of demand for final products per year'
            / (di,rf)                      .1
              (pg,rg,et,pp,bt,be,to,mx)    .15 /
   r(c,j,t) 'demand for final products   (1000 tpy)';

r(cf,j,t) = (tqcf(cf)*d(cf,j)/100.)*((1 + grfp(cf))**(midyear(t) - 1980));

display r;

Table tc(i,*) 'transportation cost   (us $ per ton)'
               inchon  seoul
   inchon                .84
   yosu          1.86   2.65
   ulsan         5.00   3.62

   +           kwangju  yosu  taegu  ulsan  pusan
   inchon         5.44  7.77   5.44   8.55   4.90
   yosu           3.54         3.00   3.72   1.77
   ulsan          5.03  3.72    .69          2.59;

* there are 4 modes of transportation: vessel, railroad tanker car (rtc), tanker truck (tt) and pipeline (pl)
* each transportation mode is tied to each demand region and the unit transportation cost for each
* transportation mode is different from each other and even for the same transportation mode, the unit
* transportation cost is different according to the total volume transported.

Table pr(c,*) 'commodity prices      (us $ per ton)'
               imports  exports
   sa            231.7
   et            908.      705.
   pp            704.      529.
   bt            706.      629.
   mx            662.      373.
   to            662.      526.
   be            706.      373.;

Scalar
   eu  'export upper bound (1000 tpy)' / 400 /
   imu 'import upper bound (1000 tpy)' / 300 /;

Parameter
   pv(c) 'import prices (us $ per ton)'
   pe(c) 'export prices (us $ per ton)';

pv(crcf) = pr(crcf,"imports");
pe(crcf) = pr(crcf,"exports");

Table inv(m,*) 'investment data'
* size   capacity limit of economy of scale (1000 tpy)
* cost   cost of the production unit at size "size" (1000 us$)
* scale  scale factor;  cost = xx*size ** scale
         size    cost  scale
   ad   10000   50000    .67
   cr    1000   40000    .67
   cc    2000  200000    .67
   ds    2000  160000    .67
   sc     500  200000    .67
   bx     100   50000    .67
   pf     200  150000    .67
   ar     400  150000    .67
   hd     200  150000    .67;

Set is 'investment function segment' / 1*4 /;

Scalar
   life  'life of productive units (years)' / 20   /
   discr 'discount rate'                    /   .1 /
   caprf 'capital recovery factor';

Parameter
   site(i)       'site factor' / inchon 1.5 , (ulsan,yosu) 1 /
   omega(m,is,i) 'cost increase for small plants (1000 us$)'
   ss(m,is)      'investment segment size        (1000 tpy)'
   discf         'discount factor';

inv(m,"fixed") = inv(m,"cost")*(.5**(inv(m,"scale") - 1) - 1);
omega(m,"1",i) = inv(m,"fixed")*site(i);
omega(m,"2",i) = inv(m,"cost")*site(i);
omega(m,"3",i) = omega(m,"2",i)*5 ;
omega(m,"4",i) = omega(m,"2",i)*10*1.25;

ss(m,"2") = inv(m,"size");
ss(m,"3") = ss(m,"2")*5;
ss(m,"4") = ss(m,"2")*10;

life     = 20;
discr    = .1;
caprf    = discr/(1 - (1 + discr)**(-life));
discf(t) = (1 + discr)**(1980 - midyear(t));

display inv, omega, ss, caprf, discf;

Table qll(c,q) 'quality bounds (lower bounds)'
        oc
   pg   95.
   rg   85.;

Table quu(c,q) 'quality bounds (upper bounds)'
        va    su
   pg   12
   rg   12
   di       1.00
   rf       4.00;

Table at 'attributes of commodities in blending'
               oc    va    su
   sa.adga   84.0  16.0
   sa.adna   59.4   2.5
   sa.adke                .26
   sa.adfo               1.02
   sa.adrs               4.35
   sa.crga   115.0  5.0
   sa.crna    97.0  2.5
   sa.ccga    93.7  6.9
   sa.ccfo                .91
   sa.dsrs               1.00;

Parameter suurf(t) 'upper limit of sulfur content in rf'
                   / 1980-84  4, 1985-89  3.5, 1990-94  3, 1995-99  2.5 /;

$sTitle Model Definition
Positive Variable
   z(c,p,i,t)    'process level                            (1000 tpy)'
   w(c,c,c,i,t)  'blending level at refinery               (1000 tpy)'
   h(m,i,t)      'capacity expansion                       (1000 tpy)'
   s(m,is,i,t)   'scale segment of investment              (1000 tpy)'
   xf(c,i,j,t)   'domestic shipment: final products        (1000 tpy)'
   xi(c,c,i,i,t) 'domestic shipment: intermediates         (1000 tpy)'
   vf(c,j,t)     'imports: final products                  (1000 tpy)'
   vr(c,i,t)     'imports: raw materials-crude oil         (1000 tpy)'
   e(c,i,t)      'exports: final products                  (1000 tpy)';

Variable
   tcost         'total cost                              (1000 us $)'
   rawmat(t)     'raw material cost                       (1000 us $)'
   operat(t)     'operating cost                          (1000 us $)'
   trans(t)      'transportation cost                     (1000 us $)'
   ccost(t)      'capital cost                            (1000 us $)'
   import(t)     'import cost                             (1000 us $)'
   export(t)     'export revenue                          (1000 us $)';

Binary Variable
   y(m,i,t)      'binary variable';

Equation
   mbr(c,i,t)    'raw material balance                     (1000 tpy)'
   mbir(c,c,i,t) 'intermediate material balance: refinery  (1000 tpy)'
   mbip(c,c,i,t) 'intermediate material balance: petrochem (1000 tpy)'
   mbfr(c,i,t)   'final material balance: refinery goods   (1000 tpy)'
   mbfp(c,i,t)   'final material balance: petrochemicals   (1000 tpy)'
   qcfl(c,q,i,t) 'quality constraints: lower bounds'
   qcfu(c,q,i,t) 'quality constraints: upper bounds'
   cc(m,i,t)     'capacity constraints                     (1000 tpy)'
   id(m,i,t)     'investment variable definition           (1000 tpy)'
   ic(m,i,t)     'combination of 0 and 1'
   mdcf(c,j,t)   'market demand of final products          (1000 tpy)'
   eub(t)        'exports constraints: upper bounds        (1000 tpy)'
   imub(t)       'imports constraints: upper bounds        (1000 tpy)'
   obj           'accounting: total cost                   (1000 us$)'
   araw(t)       'accounting: raw material cost            (1000 us$)'
   aoper(t)      'accounting: operating cost               (1000 us$)'
   atrans(t)     'accounting: transportation cost          (1000 us$)'
   acap(t)       'accounting: investment cost              (1000 us$)'
   aim(t)        'accounting: import cost                  (1000 us$)'
   aex(t)        'accounting: export revenue               (1000 us$)';

mbr(cr,i,t).. sum(p, a(cr,cr,p)*z(cr,p,i,t)) + vr(cr,i,t) =g= 0;

mbir(cr,cb,i,t)..   sum(p, a(cr,cb,p)*z(cr,p,i,t))
                 +  sum(ip$tc(i,ip), xi(cr,cb,ip,i,t) - xi(cr,cb,i,ip,t))$cs(cb)
                =g= sum(cfr$cfcb(cfr,cb), w(cr,cb,cfr,i,t));

mbip(cr,ci,i,t)..   sum(p, a(cr,ci,p)*z(cr,p,i,t))
                 +  sum(ip$tc(i,ip), xi(cr,ci,ip,i,t) - xi(cr,ci,i,ip,t))$cs(ci)
                =g= 0;

mbfr(cfr,i,t).. sum((cr,cb)$cfcb(cfr,cb), w(cr,cb,cfr,i,t)) =e= sum(j, xf(cfr,i,j,t));

mbfp(cfp,i,t).. sum((cr,p), a(cr,cfp,p)*z(cr,p,i,t)) =g= sum(j, xf(cfp,i,j,t)) + e(cfp,i,t);

qcfl(cfr,q,i,t)$qll(cfr,q)..
   sum((cr,cb)$cfcb(cfr,cb), at(cr,cb,q)*w(cr,cb,cfr,i,t)) =g= sum(j, qll(cfr,q)*xf(cfr,i,j,t));

qcfu(cfr,q,i,t)$quu(cfr,q)..
   sum((cr,cb)$cfcb(cfr,cb), at(cr,cb,q)*w(cr,cb,cfr,i,t)) =l= sum(j, quu(cfr,q)*xf(cfr,i,j,t))$cfrerf(cfr)
                                                            +  sum(j, suurf(t)*xf(cfr,i,j,t))$cfrrf(cfr);

cc(m,i,t)..    sum(p, b(m,p)*sum(cr, z(cr,p,i,t))) =l= kat1(m,i) + sum(tep$ts(t,tep), h(m,i,tep));

id(m,i,te)..   h(m,i,te) =e= sum(is, ss(m,is)*s(m,is,i,te));

ic(m,i,te)..   y(m,i,te) =e= sum(is, s(m,is,i,te));

mdcf(cf,j,t).. sum(i, xf(cf,i,j,t)) + vf(cf,j,t)$cfp(cf) =g= r(cf,j,t);

eub(t)..       sum((cfp,i), e(cfp,i,t))  =l= eu;

imub(t)..      sum((cfp,j), vf(cfp,j,t)) =l= imu;

obj..          tcost =e= year*sum(t, discf(t)*(rawmat(t) + operat(t) + trans(t) + ccost(t) + import(t) - export(t)));

araw(t)..      rawmat(t) =e= sum((cr,i), pv(cr)*vr(cr,i,t));

aoper(t)..     operat(t) =e= sum(p, op(p)*sum((cr,i), z(cr,p,i,t)));

atrans(t)..    trans(t)  =e= sum((cf,i,j), tc(i,j)*xf(cf,i,j,t)) + sum((cr,cs,i,ip), tc(i,ip)*xi(cr,cs,i,ip,t));

acap(t)..      ccost(t)  =e= caprf*sum((tep,m,is,i)$ts(t,tep), omega(m,is,i)*s(m,is,i,tep));

aim(t)..       import(t) =e= sum((cfp,j), pv(cfp)*vf(cfp,j,t));

aex(t)..       export(t) =e= sum((cfp,i), pe(cfp)*e(cfp,i,t));

Model petro / all /;

solve petro minimizing tcost using mip;

$sTitle Report
Parameter
   vrlev(i,t)     'crude oil imports level            (1000 ton)'
   wlev(c,c,i,t)  'blending level: detailed           (1000 ton)'
   wlevf(c,i,t)   'blending level: aggreated          (1000 ton)'
   tzlev(p,i,t)   'total production level             (1000 ton)'
   tzplev(p,t)    'total production per process       (1000 ton)'
   xilev(c,i,i,t) 'shipment of intermediate products  (1000 ton)'
   xflev(c,j,t)   'shipment of final products         (1000 ton)'
   tic(m,i,t)     'total installed capacity per plant (1000 ton)'
   tict(m,t)      'total installed capacity           (1000 ton)'
   tac(m,t)       'total added capacity               (1000 ton)'
   ti(m,i,t)      'total investment cost per unit    (1000 us $)'
   tit(i,t)       'total investment cost per plant   (1000 us $)'
   tcu            'undiscounted total cost           (1000 us $)';

vrlev(i,t)       = sum(cr, vr.l(cr,i,t));
wlev(cb,cf,i,t)  = sum(cr, w.l(cr,cb,cf,i,t));
wlevf(cf,i,t)    = sum((cr,cb), w.l(cr,cb,cf,i,t));
tzlev(p,i,t)     = sum(cr, z.l(cr,p,i,t));
tzplev(p,t)      = sum((cr,i), z.l(cr,p,i,t));
xilev(ci,ip,i,t) = sum(cr, xi.l(cr,ci,ip,i,t));
xflev(cf,j,t)    = sum(i, xf.l(cf,i,j,t));
tic(m,i,t)       = kat1(m,i) + sum(tp$ts(t,tp), h.l(m,i,tp));

tict(m,t) = sum(i, tic(m,i,t));
tac(m,t)  = sum(i, h.l(m,i,t));
ti(m,i,t) = sum(is, omega(m,is,i)*s.l(m,is,i,t));
tit(i,t)  = sum(m, ti(m,i,t));
tcu       = sum(t, year*(rawmat.l(t) + operat.l(t) + trans.l(t) + ccost.l(t) + import.l(t) - export.l(t)));

display vrlev, wlev, wlevf, tzlev, tzplev, xilev;
display xf.l, xflev, vf.l, e.l, h.l, s.l, tac, tic, tict, ti, tit, tcu;