!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ
!   A two-stage stochastic mixed integer linear program for optimal bidding of hydropower into Nord Pool Spot                 İİİİİİİİİ
!   Developed by Erik Alnĉs and Roger Grĝndahl, NTNU                                                                          İİİİİİİİİ
!   Last updated June 7, 2012                                                                                                 İİİİİİİİİ
!   Constraints are numbered according to mathematical formulation in "Bidding Revealed - An empirical analysis"
!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ

model BidOptimizer
uses "mmxprs";                                                                                   !Access to the Xpress-Optimizer solver
uses "mmsystem";                                                                                    !Access to operating system actions
uses "mmive";                                                                                !Access to visual environment for graphing
uses "mmodbc";
!uses "mmquad";                                                                                      !Access to external ODBC databases

options explterm                                                                                       !Expressions always end with ";"
options noimplicit                                                                                       !All declarations are explicit

parameters
    Input   =   "mmodbc.excel:noindex;Input_cascadeXX.xlsx";
    Results =   "mmodbc.excel:noindex;Results.xlsx";
    Results2=   "mmodbc.excel:noindex;Input_cascade(XX+1).xlsx";
end-parameters

declarations
    start_time:     real;
    end_time:       real;
    running_time:   real;
end-declarations

start_time:=gettime;
setparam("XPRS_MIPRELSTOP",0.02);
setparam("XPRS_VERBOSE",TRUE);
!setparam("XPRS_PRESOLVE",0);  !Through tuning we found that turning presolve off might decrease runtime drastically for some instances

!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ
!   Sets                                                                                                                      İİİİİİİİİ
!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ

declarations
!Day-ahead
nH,nHd,nS,nI,nB:                            integer;
S:                  set                 of  integer;                                                                  !Set of scenarios
H_d:                set                 of  integer;                                                        !Hours of day-ahead bidding
I:                  set                 of  integer;                                                                  !Set of bidpoints
B:                  set                 of  integer;                                                     !Set of possible blocks to bid
R:                  set                 of   string;                                                                        !Reservoirs
T:                  set                 of   string;                                                                          !Turbines
K:                  set                 of   string;                                                                    !Power Stations
!Rest of week
H:                  set                 of  integer;                                                    !All hours the model spans over
H_l:                set                 of  integer;                                              !The hours that are run deterministic
end-declarations

initializations from Input
nHd         as          "nHours";
nS          as          "nScenarios";
nI          as          "nBidpoints";
nB          as          "nBlocks";
R           as          "Reservoirs";
T           as          "Turbines";
K           as          "PowerStations";
!Rest of week
nH          as          "nHoursLong";
end-initializations

H_d     :=  1 .. nHd;
H_l     :=  nHd+1 .. nH;
H       :=  1 .. nH;
S       :=  1 .. nS;
I       :=  1 .. nI;
B       :=  1 .. nB;

finalize(H_d);
finalize(H_l);
finalize(H);
finalize(S);
finalize(I);
finalize(B);

!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ
!   Coefficients                                                                                                              İİİİİİİİİ
!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ
declarations
!Day-ahead
T_in_K:             array(K) of set of string; !All turbines in power station K
K_in_T:             array(T)        of string;
K_in_R:             array(K)        of string;
K_above_R:          array(K)        of string;
nE:                 array(T)        of integer;
E:                  array(T)        of set of integer;
W:                  dynamic array(T,1..6)   of real;
W0:                 dynamic array(T,1..6)   of real;

Maxwatt:            array(T,H)      of real;
MinWatt:            array(T)        of real;
MinFlow:            array(K,H)      of real;
SPOT:               array(S,H)      of real;
F:                  array(R,H)      of real;                                                                                    !Inflow
PROB:               array(S)        of real;                                                                             !Probabilities
P:                  array(I)        of real;                                                                       !Price bid for bid I 
C:                  array(T)        of real;                                                                            !Start-up costs
C_feedVar:          array(K,H)      of real;                                                          !Variable part of the feed-in fee
C_fixed:            real;                                                                                !Fixed part of the feed-in fee
ResStations:        array(R)        of set of string;                                 !Power stations that draws water from reservoir R
ResInitial:         array(R)        of real;
ResEnd:             array(R)        of real;
ResMax:             array(R)        of real;
ResMin:             array(R)        of real;
Potential:          array(T)        of real;
WaterUsed:          array(K)        of real;

temp:               string;
Time:               array(R,K)      of integer;                                 !Time for water to go from reservoir to a power station
StationsAbove:      array(R)        of set of string;                                   !Power stations from which water flows to res R

BlockStart:         array(B)        of integer;
BlockEnd:           array(B)        of integer;
P_b:                array(B,I)      of real;
BlockLength:        array(B)        of integer;
BlockAverage:       array(S,B)      of real;
end-declarations

initializations from Input
K_in_T;
K_in_R;
K_above_R;
W;
W0;
nE;
SPOT                as          "Prices";
F                   as          "Inflow";
PROB                as          "Probabilities";
P                   as          "BidPrices";
Maxwatt             as          "MaxWatt";
MinWatt;
MinFlow             as          "MinFlow";
C                   as          "Startcost";
C_feedVar           as          "C_feedIn";
C_fixed;
ResInitial          as          "ResInitial";
ResEnd;
ResMin;             
Potential;
ResMax              as          "ResMax";
BlockStart;
BlockEnd;
WaterUsed;
end-initializations

forall(s in S, b in B)do
    BlockLength(b):=0;
    BlockAverage(s,b):=0;
        forall( h in H| h>=BlockStart(b) and h<=BlockEnd(b))do
        BlockAverage(s,b)   :=BlockAverage(s,b) + SPOT(s,h);
        BlockLength(b)      :=BlockLength(b) + 1;
        end-do
    BlockAverage(s,b) := BlockAverage(s,b) / BlockLength(b);
end-do

forall(b in B,i in I)do
P_b(b,i):=P(i);
end-do

forall(t in T)do
E(t):=1..nE(t); 
end-do

forall(t in T)do
temp:=K_in_T(t);
T_in_K(temp)+={t};
end-do

forall(k in K)do
temp:=K_in_R(k);
ResStations(temp)+={k};
end-do

forall(k in K|K_above_R(k)<>"")do
temp:=K_above_R(k);
StationsAbove(temp)+={k};
end-do

forall(r in R) do
    forall(k in StationsAbove(r))do
    Time(r,k):=0;                                      !Set to zero if the travel time between reservoirs are shorter than half an hour
    end-do
end-do

!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ
!   Decision variables                                                                                                        İİİİİİİİİ
!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ

declarations
x_h:        array(I,H_d)            of mpvar;                                                           !Volume bid for bid I in hour H
y_h:        array(S,H_d)            of mpvar;                                      !Commitment for hourly bids in scenario S and hour H
x_b:        array(B,I)              of mpvar;                                                               !Volume bid for block bid B
x_bh:       array(H)                of mpvar;                                                       !Volume bid as block bids in hour H
y_b:        array(S,B)              of mpvar;                                       !Commitment for block bids in scenario S and hour H
y:          array(S,H_d)            of mpvar;                                                                   !Total power commitment
watt:       array(S,H,T)            of mpvar;                                                              !Power output from turbine T
q:          array(S,H,T)            of mpvar;                                                                 !Water usage in turbine T
q_k:        array(S,H,K)            of mpvar;                                                           !Water usage in power station K
running:    array(S,H_d,T)          of mpvar;                              !Binary variable stating whether turbine T is running or not
started:    array(S,H_d,T)          of mpvar;                       !Binary variable stating whether turbine T started or not in hour H
resLevel:   array(S,H,R)            of mpvar;                                                    !Level of MWh in reservoir at hour end
end-declarations

forall(s in S, h in H_d, t in T) do
running(s,h,t) is_binary;
started(s,h,t) is_binary;
end-do

!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ
!   Constraints                                                                                                               İİİİİİİİİ
!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ

declarations
IncreasingCurves:   array(S,H_d,I)          of linctr;                                                                         !Eq.(4.4)
BidCommitment:      array(S,H_d,I)          of linctr;                                                                         !Eq.(4.5)
BlockCommitment:    array(S,B)              of linctr;                                                                         !Eq.(4.6)
BlockBidHour:       array(H_d)              of linctr;                                                                        
BidHour:            array(H_d)              of linctr;                                                                         !Eq.(4.7)
YSum:               array(S,H_d)            of linctr;                                                                         !Eq.(4.8)
WattSum:            array(S,H_d)            of linctr;                                                                         !Eq.(4.9)
RunningCon:         array(S,H_d,T)          of linctr;                                                                        !Eq.(4.10)
RunningCon2:        array(S,H_d,T)          of linctr;                                                                        !Eq.(4.11)
StartedCon:         array(S,H_d,T)          of linctr;                                                                        !Eq.(4.12)
LoopCon:            array(S,T)              of linctr;                                                                        !Eq.(4.13)
MaxCap:             array(S,H,T)            of linctr;                                                                        !Eq.(4.14)
WaterUsage:         array(S,H,T,range)      of linctr;                                                                        !Eq.(4.15)
QSum:               array(S,H,K)            of linctr;                                                                        !Eq.(4.16)
MinFlowCon:         array(S,H,K)            of linctr;                                                                        !Eq.(4.17)
ResBalance:         array(S,H,R)            of linctr;                                                                 !Eq.(4.18)-(4.19)
ResMaxBalance:      array(S,H,R)            of linctr;                                                                        !Eq.(4.20)
ResMinBalance:      array(S,H,R)            of linctr;                                                                        !Eq.(4.21)

!WaterCoupling:     array(S,R)              of linctr;                                                                        !Eq.(4.22)
WaterCoupling:      array(S,K)              of linctr;                                                                        !Eq.(4.32)

end-declarations

forall(s in S, h in H_d, i in I|i<>1)do
IncreasingCurves(s,h,i):=
x_h(i-1,h) <= x_h(i,h);
end-do

forall(s in S, h in H_d, i in I|(i<>1 and SPOT(s,h)>P(i-1) and SPOT(s,h)<P(i))) do
BidCommitment(s,h,i):=
y_h(s,h) = x_h(i-1,h)+ (SPOT(s,h)-P(i-1))*(x_h(i,h)-x_h(i-1,h))/((P(i)-P(i-1)));                    
end-do

forall(s in S, b in B)do
BlockCommitment(s,b):=
y_b(s,b)= sum(i in I|P_b(b,i)<=BlockAverage(s,b)) x_b(b,i);
end-do

forall(h in H_d)do
BlockBidHour(h):=
x_bh(h)=sum(i in I,b in B|BlockStart(b)<=h and BlockEnd(b)>=h) x_b(b,i);
end-do

forall(h in H_d)do  
BidHour(h):=
x_h(nI,h)+x_bh(h)<=sum(t in T)Maxwatt(t,h);
end-do

forall(s in S,h in H_d)do
YSum(s,h):=
y(s,h) = y_h(s,h) + sum(b in B|BlockStart(b)<=h and BlockEnd(b)>=h) y_b(s,b); 
end-do

forall(s in S,h in H_d)do
WattSum(s,h):=
sum(t in T)watt(s,h,t) = y(s,h);
end-do

forall(s in S, h in H_d, t in T) do
RunningCon(s,h,t):=
watt(s,h,t)<=running(s,h,t)*Maxwatt(t,h);
end-do

forall(s in S, h in H_d, t in T) do
RunningCon2(s,h,t):=
watt(s,h,t)>=running(s,h,t)*MinWatt(t);
end-do

forall(s in S, h in H_d, t in T|h<>1) do
StartedCon(s,h,t):=
started(s,h,t) >= running(s,h,t) - running(s,h-1,t);
end-do

forall (s in S, t in T) do     !This could alternatively be replaced by a discount rate so that 
LoopCon(s,t):=                 !start-ups in hour 23 cost less than in hour 2   
started(s,1,t) >= running(s,1,t) - running(s,nHd,t);
end-do

forall(s in S, h in H_l,t in T)do
MaxCap(s,h,t):=
watt(s,h,t)<=Maxwatt(t,h);
end-do

forall(s in S, h in H, t in T, e in E(t)) do
WaterUsage(s,h,t,e):=
watt(s,h,t)<=W0(t,e) + W(t,e)*q(s,h,t);
end-do

forall(s in S, h in H, k in K)do
QSum(s,h,k):=
q_k(s,h,k)=sum(t in T_in_K(k))q(s,h,t);
end-do

forall(s in S, h in H, k in K)do
MinFlowCon(s,h,k):=
q_k(s,h,k)>=MinFlow(k,h);
end-do

forall(s in S,r in R)do
ResBalance(s,1,r):=
resLevel(s,1,r)=ResInitial(r) - sum(k in ResStations(r))q_k(s,1,k) + F(r,1);
end-do

forall(s in S,r in R, h in H|h<>1)do
ResBalance(s,h,r):=
resLevel(s,h,r)=resLevel(s,h-1,r)   - sum(k in ResStations(r))q_k(s,h,k) + F(r,h) + sum(k in StationsAbove(r)) q_k(s,h-Time(r,k),k);
end-do

forall(s in S,r in R, h in H)do
ResMaxBalance(s,h,r):=
resLevel(s,h,r)<= ResMax(r);
end-do

forall(s in S,r in R, h in H)do
ResMinBalance(s,h,r):=
resLevel(s,h,r)>= ResMin(r);
end-do


!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ
!   Deterministic long term part                                                                                              İİİİİİİİİ
!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ
(!forall(s in S,r in R)do
WaterCoupling(s,r):=
resLevel(s,nH,r)=ResEnd(r);                                                            !This couples the reservoir levels at period end
end-do!)

forall(s in S, k in K)do
WaterCoupling(s,k):=                                                                !This can be used instead of the reservoir coupling 
sum(h in H)q_k(s,h,k)<=WaterUsed(k);
end-do



!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ
!   Objective                                                                                                                 İİİİİİİİİ
!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ

declarations
SpotIncome:     array(H_d) of linctr;
StartCost:      array(H_d) of linctr;
FeedInCost:     array(H_d) of linctr;
DailyProfit:    linctr;

!Long term
SpotProfit_long:                    linctr;
FeedInCost_long:                    linctr;
LongProfit:                         linctr;
end-declarations

forall(h in H_d)do
SpotIncome(h):=     sum(s in S) y(s,h)*SPOT(s,h)*PROB(s);
StartCost(h):=      sum(s in S, t in T) started(s,h,t)*C(t)*PROB(s);
FeedInCost(h):=     sum(s in S, k in K, t in T_in_K(k)) watt(s,h,t)*(C_fixed+C_feedVar(k,h)*SPOT(s,h))*PROB(s);
end-do

DailyProfit:=       sum(h in H_d)( SpotIncome(h) - StartCost(h)- FeedInCost(h) );

!long term
SpotProfit_long:=   sum(s in S, h in H_l, k in K, t in T_in_K(k))( PROB(s)*watt(s,h,t)*(SPOT(s,h)*(1-C_feedVar(k,h))-C_fixed));
LongProfit:= DailyProfit + 1.007*SpotProfit_long;                                                                     !This is the Beta

maximize(LongProfit);
end_time:=gettime;
running_time:=end_time-start_time;



!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ
!   Postprocessing                                                                                                            İİİİİİİİİ
!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ

declarations
Print_x:            array(I,H_d)                of real;
Print_spot:         array(H_d)                  of real;
Print_cost:         array(H_d)                  of real;
Print_Feed:         array(H_d)                  of real;
Print_q:            array(H_d)                  of real;
Print_pot:          array(H_d)                  of real;
Print_watt:         array(S,H)                  of real;
Print_longObjective:                               real;
Print_block:        dynamic array(B)            of real;
Print_block_price:  dynamic array(B,I)          of real;
Print_block_start:  dynamic array(B)            of real;
Print_block_end:    dynamic array(B)            of real;
Print_block_volume: dynamic array(B,I)          of real;
Print_block_com:    dynamic array(B)            of real;
Print_Reservoir:    array(R)                    of real;
Print_waterUsed:    real;
Print_priceForecast:array(H_d)                  of real;
waterOutput:        array(K)                    of real;
end-declarations

forall(i in I, h in H_d)do
Print_x(i,h):=getsol(x_h(i,h)); 
end-do

forall(h in H_d)do
Print_spot(h)           :=getsol(SpotIncome(h));
Print_cost(h)           :=getsol(StartCost(h));
Print_Feed(h)           :=getsol(FeedInCost(h));
Print_pot(h)            :=sum(s in S, t in T) getsol(q(s,h,t)*Potential(t)*PROB(s));
Print_q(h)              :=sum(s in S, k in K)getsol(q_k(s,h,k))*PROB(s);
Print_priceForecast(h)  :=SPOT(1,h);
end-do

forall(s in S,h in H)do
Print_watt(s,h) :=sum(t in T)getsol(watt(s,h,t));
end-do

Print_longObjective:=getsol(SpotProfit_long);

forall(b in B, i in I)do
    if(getsol(x_b(b,i))<>0)then
    Print_block(b)          :=b;
    Print_block_price(b,i)  :=P_b(b,i);
    Print_block_start(b)    :=BlockStart(b);
    Print_block_end(b)      :=BlockEnd(b);
    Print_block_volume(b,i) :=getsol(x_b(b,i));
    Print_block_com(b)      :=getsol( sum(s in S)( y_b(s,b)*PROB(s)) );
    end-if
end-do

forall(r in R)do
Print_Reservoir(r):=sum(s in S) PROB(s)*getsol(resLevel(s,nHd,r));
end-do

Print_waterUsed := sum(s in S, k in K, h in H_d)getsol(q_k(s,h,k))*PROB(s);

forall(k in K)do
waterOutput(k) := sum(s in S, h in H_d)getsol(q_k(s,h,k))*PROB(s);
end-do

initializations to Results
nS;
nH;
Print_spot         as  "SpotIncome";
Print_cost         as  "StartCost";
Print_Feed         as  "FeedCost";
Print_watt         as  "Watt";
Print_block        as  "B";
Print_block_price  as  "BlockPrice";
Print_block_start  as  "BlockStart";
Print_block_end    as  "BlockEnd";
Print_block_volume as  "BlockVolume";
Print_block_com    as  "BlockCom";
Print_longObjective;
running_time;
P;
Print_x            as  "BidCurve";
Print_q;
Print_pot;
Print_waterUsed;
Print_priceForecast;
end-initializations

initializations to Results2
Print_Reservoir    as  "ResInitial";
waterOutput;
end-initializations

!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ
end-model                                                                                                                    !İİİİİİİİİ
!İİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİİ