* MARKET CLEARING UNDER EQUIPMENT FAILURES

Sets
N nodes     /n1*n3/
I units     /i1*i3/
J loads     /j1/
T periods   /t1*t4/
w scenarios /w0*w8/
M offer blocks  /m1/
S system states /s0*s2/;

Alias(n,r);

scalars
nt number of time periods /4/
dt duration of time periods /1/
x  reactance of transmission lines /0.13/
Pbase  power base /41/;

Table gdata(I,*) generation data
      Pmin     pmax      Istatus    pini
i1    10       100       0          0
i2    10       100       0          0
i3    10       50        0          0   ;

Table rdata_g(i,t,*) generation-side reserve data
             rspin_max_up        rspin_max_dn       rnspin_max
i1.t1*t4     90                  90                 100
i2.t1*t4     90                  90                 100
i3.t1*t4     40                  40                 50    ;


Table    costblock(i,t,m) offer cost of energy blocks
              m1
i1.t1*t4      30
i2.t1*t4      40
i3.t1*t4      20 ;

Table widthblock(i,t,m)   width of energy blocks
              m1
i1.t1*t4      100
i2.t1*t4      100
i3.t1*t4      50 ;

Table cr_schg(i,t,*)   generation-side reserve offer cost
               up      dn     ns
i1.t1*t4       5       5      4.5
i2.t1*t4       7       7      5.5
i3.t1*t4       8       8      7        ;
 table ml(j,n) mapping of the set of loads into the set of nodes
         n1        n2      n3
j1       0         0       1    ;

Table netstate(n,n,s)  network states
        s0  s1  s2
n1.n1   0   0   0
n1.n2   1   1   1
n1.n3   1   1   0
n2.n1   1   1   1
n2.n2   0   0   0
n2.n3   1   1   1
n3.n1   1   1   0
n3.n2   1   1   1
n3.n3   0   0   0   ;

Table genstate(i,n,s)  states of units avaibility
        s0  s1  s2
i1.n1   1   1   1
i1.n2   0   0   0
i1.n3   0   0   0
i2.n1   0   0   0
i2.n2   1   1   1
i2.n3   0   0   0
i3.n1   0   0   0
i3.n2   0   0   0
i3.n3   1   0   1   ;

parameters
elast_limits(j,t,*)  elasticity limits of loads
/j1.t1.sup     30
 j1.t2.sup     80
 j1.t3.sup     110
 j1.t4.sup     40
 j1.t1.inf     30
 j1.t2.inf     80
 j1.t3.inf     110
 j1.t4.inf     40
/

Rdata_d(j,t,*)   demand-side reserve data
/j1.t1.rspin_max_up     3
 j1.t2.rspin_max_up     8
 j1.t3.rspin_max_up     11
 j1.t4.rspin_max_up     4
 j1.t1.rspin_max_dn     3
 j1.t2.rspin_max_dn     8
 j1.t3.rspin_max_dn     11
 j1.t4.rspin_max_dn     4
/

prob(w)  scenario w probability
/w0       0.98807312
 w1       0.00199015
 w2       0.00199015
 w3       0.00199015
 w4       0.00199015
 w5       0.00099157
 w6       0.00099157
 w7       0.00099157
 w8       0.00099157
/

TAU(w)   period in which the contigency defining scenario w occures
/w0       5
 w1       1
 w2       2
 w3       3
 w4       4
 w5       1
 w6       2
 w7       3
 w8       4
/


num0(i,t)         number of energy blocks offered by each unit
lambdaSU(i,t)      start-up offer cost
lambdal(j,t)       demand utility
cr_schd(j,t,*)    demand-side reserve offer cost
vlol(j,t)         value of lost load
b(n,n)            imaginary parts of the admittance of lines
fmax(n,n)         transmission capacity limits
conec(n,n,t,w)    in time period t and scenario w
MG(i,n,t,w)       mapping of the set of units into the set of nodes
G(i,t,w)          set of available units per scenario and period ;

num0(i,t)=1;
lambdasu(i,t)=100;
lambdal(j,t)=0;
b(n,r)$(ord(n) ne ord(r))= -1/x ;
G(i,t,w)=1;
cr_schd(j,t,'up')=70;
cr_schd(j,t,'dn')=70;
vlol(j,t)=1000;
fmax(n,r)$(ord(n) ne ord(r))= 55;

loop(t,
     loop(w,
          if(ord(T)<TAU(w),
              conec(n,r,t,w)=netstate(n,r,'s0');
              MG(i,n,t,w)=genstate(i,n,'s0');
          else
             if(ord(w)<=5,
                conec(n,r,t,w)=netstate(n,r,'s1');
                MG(i,n,t,w)=genstate(i,n,'s1');
                G('i3',t,w)=0;
               else
                conec(n,r,t,w)=netstate(n,r,'s2');
                mg(i,n,t,w)=genstate(i,n,'s2');
            );
                );
                    );
                        );
Variables

EC            expected cost(objective function value)

*First-stage variables

Ps(i,t)                power output schedule
Ls(j,t)                demand schedule
r_spin_schg_d(i,t)     generation-side spinning reserve schedule (downward)
r_spin_schg_u(i,t)     generation-side spinning reserve schedule (upward)
r_spin_schd_d(j,t)     demand-side spinning reserve schedule (downward)
r_spin_schd_u(j,t)     demand-side spinning reserve schedule (upward)
r_nspin_schg(i,t)      generation-side non-spinning reserve schedule
u(i,t)                 scheduled commitment status

*Second-stge varibles

Lc(j,t,w)              power consumption in real time
Lshed(j,t,w)           amount of involuntarily load shed
Pg(i,t,w)              generator power output in real time
pgblock(i,t,w,m)       power produced from energy blocks
r_spin_depg_d(i,t,w)   generation-side spinning reserve deployment (downward)
r_spin_depg_u(i,t,w)   generation-side spinning reserve deployment (upward)
r_spin_depd_d(j,t,w)   demand-side spinning reserve deployment (downward)
r_spin_depd_u(j,t,w)   demand-side spinning reserve deployment (upward)
r_nspin_depg(i,t,w)    generation-side non-spinning reserve deployment
angle(n,t,w)           vlotage angle
pinj(n,t,w)            power injection
f(n,r,t,w)             power flow
v(i,t,w)               real-time commitment status
Csu(i,t,w)             start-up cost ;

positive variables
Csu(i,t,w), Lshed(j,t,w), r_spin_depg_d(i,t,w), r_spin_depg_u(i,t,w),r_spin_depd_d(j,t,w),  r_spin_depd_u(j,t,w),r_nspin_depg(i,t,w),pgblock(i,t,w,m), r_nspin_schg(i,t), r_spin_schg_d(i,t) , r_spin_schg_u(i,t) ,r_spin_schd_u(j,t),r_spin_schd_d(j,t);

Binary variables   v(i,t,w),u(i,t) ;

* Elacticity limits od demand
Ls.up(j,t)=elast_limits(j,t,'sup');
Ls.lo(j,t)=elast_limits(j,t,'inf');

angle.fx('n1',t,w)=0;

lshed.fx(j,t,w)$(ord(T) LT TAU(w))=0;
r_spin_depg_d.fx(i,t,w)$(ord(T) LT TAU(w))=0;
r_spin_depg_u.fx(i,t,w)$(ord(T) LT TAU(w))=0;
r_spin_depd_d.fx(j,t,w)$(ord(T) LT TAU(w))=0;
r_spin_depd_u.fx(j,t,w)$(ord(T) LT TAU(w))=0;
r_nspin_depg.fx(i,t,w)$(ord(T) LT TAU(w))=0;

Equations

ECfunction      objective function

**Electricity market constraints

* production limits
Max_prod(i,t)
Min_prod(i,t)

*market equilibria
Market_equ(t)

** scheduled  reserv deterination constraints
* generation side
reserve_sch_spin_g_u(i,t)      spinning up
reserve_sch_spin_g_d(i,t)      spinning down
reserve_sch_nspin_g(i,t)       non-spinning
**demand side
reserve_sch_spin_d_u(j,t)      spinning up
reserve_sch_spin_d_d(j,t)      spinning down

**Real-time operating cnstraints
SUo(i,w)                       T=1
SU(i,t,w)                      T>1

*generation limits

GL1(i,t,w)
GL2(i,t,w)
GL3(i,t,w,M)
GL4(i,t,w)

*involuntary load shedding contraints
Limit_ENS(J,T,W)

*network constraints
PB(n,t,w)                    power balance
Def_pinj(n,t,w)              definition of power injection
Def_flow(n,r,t,w)            definition of power flows
Max_cap(n,r,t,w)             transmission capacity

*linking constraints
pgenerated(i,t,w)            composition of generator power outputs
pconsumed(j,t,w)              composition of the power consumption

*generation side reserve
reserve_dpl_spin_G_u(i,t,w)  spinning up
reserve_dpl_spin_G_d(i,t,w)  spinning down
reserve_dpl_nspin_G(i,t,w)   non-spinning
*demand side reserve
reserve_dpl_spin_d_u(j,t,w)  spinning up
reserve_dpl_spin_d_d(j,t,w)  spinning down ;

ECfunction..            EC=e=dt*sum((i,t),cr_schg(i,t,'up')*r_spin_schg_u(i,t)+cr_schG(i,t,'dn')*r_spin_schg_d(i,t)+cr_schg(i,t,'ns')*r_nspin_schg(i,t))+dt*sum((j,t),cr_schd(j,t,'up')*r_spin_schd_u(j,t)+cr_schd(j,t,'dn')*r_spin_schd_d(j,t))+dt*sum((t,w),prob(w)*(sum((i,m)$(ord(m) LE num0(i,t)), costblock(i,t,m)*pgblock(i,t,w,m))-sum(j,lambdal(j,t)*Lc(j,t,w))+sum(j,vlol(j,t)*Lshed(j,t,w))))+sum((i,t,w),prob(w)*csu(i,t,w));

Min_prod(i,t)..         ps(i,t)=g= gdata(i,'pmin')*u(i,t);
max_prod(i,t)..         ps(i,t)=l= gdata(i,'pmax')*u(i,t);
market_equ(t)..         sum(i,ps(i,t))=e=sum(j,ls(j,t));

reserve_sch_spin_g_u(i,t)..     r_spin_schg_u(i,t)=l=rdata_g(i,t,'rspin_max_up')*u(i,t);
reserve_sch_spin_g_d(i,t)..     r_spin_schg_d(i,t)=l=rdata_g(i,t,'rspin_max_dn')*u(i,t);
reserve_sch_nspin_g(i,t)..      r_nspin_schg(i,t)=l=rdata_g(i,t,'rnspin_max')*(1-u(i,t));

reserve_sch_spin_d_u(j,t)..     r_spin_schd_u(j,t)=l=rdata_d(j,t,'rspin_max_up');
reserve_sch_spin_d_d(j,t)..     r_spin_schd_d(j,t)=l=rdata_d(j,t,'rspin_max_dn');

SUo(i,w)..                      csu(i,'t1',w)=g=lambdasu(i,'t1')*(v(i,'t1',w)-gdata(i,'Istatus'));
SU(i,t,w)$(ord(t) gt 1)..       csu(i,t,w)=g= lambdasu(i,t)*(v(i,t,w)-v(i,t-1,w));

GL1(i,t,w)..                    pg(i,t,w)=g= gdata(i,'pmin')*v(i,t,w);
GL2(i,t,w)..                    pg(i,t,w)=g= gdata(i,'pmax')*v(i,t,w);
GL3(i,t,w,m)$(ord(m) LE num0(i,t))..     pgblock(i,t,w,m)=l=widthblock(i,t,m) ;
GL4(i,t,w)..                   pg(i,t,w)=e=sum(M$(ord(m) LE num0(i,t)),pgblock(i,t,w,m));

Limit_ens(j,t,w)..              lshed(j,t,w)=l= lc(j,t,w);

PB(n,t,w)..                     sum(i$(MG(i,n,t,w) eq 1), pg(i,t,w))-sum(j$(ML(j,n) eq 1), Lc(j,t,w)-lshed(j,t,w))-pinj(n,t,w)=e=0 ;

Def_pinj(n,t,w)..               pinj(n,t,w)=e= sum(r$(conec(n,r,t,w) eq 1), f(n,r,t,w));

Def_Flow(n,r,t,w)$(conec(n,r,t,w) eq 1)..   f(n,r,t,w)=e= -b(n,r)*(angle(n,t,w)-angle(r,t,w));

Max_cap(n,r,t,w)$(conec(n,r,t,w) eq  1)..   f(n,r,t,w)=l=fmax(n,r);

pgenerated(i,t,w)$(g(i,t,w) eq 1)..         pg(i,t,w)=e= ps(i,t)+r_spin_depg_u(i,t,w)+r_nspin_depg(i,t,w)-r_spin_depg_d(i,t,w);
pconsumed(j,t,w)..                          Lc(j,t,w)=e= ls(j,t)-r_spin_depd_u(j,t,w)+r_spin_depd_d(j,t,w);

reserve_dpl_spin_g_u(i,t,w)..      r_spin_depG_u(i,t,w)=l=r_spin_schg_u(i,t);
reserve_dpl_spin_g_d(i,t,w)..      r_spin_depG_d(i,t,w)=l=r_spin_schg_d(i,t);

reserve_dpl_nspin_g(i,t,w)..       r_nspin_depG(i,t,w)=l=r_nspin_schg(i,t);

reserve_dpl_spin_d_u(j,t,w)..      r_spin_depd_u(j,t,w)=l=r_spin_schd_u(j,t);
reserve_dpl_spin_d_d(j,t,w)..      r_spin_depd_d(j,t,w)=l=r_spin_schd_d(j,t);

Model  MC   /all/;

MC.optcr=0;
option   iterlim=1e8;
option reslim= 1e10 ;

option mip=cplex;

solve MC using MIP minimizing EC;

display ec.l, ps.l;

*$libinclude matout  ec.l

*$libinclude matout  ps.l   i  t

