!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!    Operational Tidal Lagoon model                          !!!!    
!!!!    (c) 2016 Trine Rollefsen Nęss and Linn Emelie Schäffer  !!!!
!!!!    author: T. R. Nęss, rev. June 2016                      !!!!
!!!!    author: L. E. Schäffer, rev June 2016                   !!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
model Operational
    
uses "mmxprs"; !gain access to the Xpress-Optimizer solver
uses "mmsystem"; !text to int
options explterm, noimplicit;

parameters
    !input values
    PriceCase = 1;  !Price cases: 1 (current price char) or 2 (high-vol prices)
    Markets = 2;    !Number of markets included: 1 or 2
    nScenarios = 15;
    nTurbines = 2;
    Parts = 40;     
    nCycles = 1409;
    nFlowBreakPoints = 4;
    nHeadBreakPoints = 5;
    MaxTime = 4*3600;   ! Max runtime in seconds
    ScenID = 2;             !scenario tree ID for stability testing
    
    
    Prob=1/nScenarios;
    nameapp='';
    
    ProblemFile = "results\\ProblemOverview_Case" + PriceCase + "_" + Markets + ".txt"; !Output file 
    ProfitFile = "results\\Profit_Case" + PriceCase + "_" + Markets + ".txt";
    PowerFile = "results\\Power_Case" + PriceCase + "_" + Markets + ".txt";
    SalesFile = "results\\Sales_Case" + PriceCase + "_" + Markets + ".txt";
    BidFile = "results\\Bids_Case" + PriceCase + "_" + Markets + ".txt";
    
    DataFolder = "data\\";
    
    Data_file = DataFolder + 'Data.txt';
    Cycle_file = DataFolder + 'C' + nCycles + '_FP' + nFlowBreakPoints + '_HP' +nHeadBreakPoints + '.txt';
    Price_file = DataFolder + 'PriceS' + nScenarios + '_' + ScenID + '.txt';
    RedCycle_file = DataFolder + 'RedCycleData' + nameapp + '.txt';
    Power_file = DataFolder + 'MaxPower.txt';
end-parameters

if(MaxTime>0.0) then
    setparam("XPRS_maxtime", MaxTime);
end-if

!Algorithm control parameters
setparam("XPRS_cutstrategy", 3);!Cutting strategy. Values:
                                    !-1: Automatic selection (default)
                                    ! 0: No cuts
                                    ! 1: Conservative cut strategy
                                    ! 2: Moderate cut strategy
                                    ! 3: Aggressive cut strategy
                                    
setparam("XPRS_heurstrategy", 3);!Heuristic strategy
                                    !-1
                                    ! 0
                                    ! 1
                                    ! 2
                                    ! 3: Extensive heuristic strategy
                                                                        

declarations
    status:array({XPRS_OPT,XPRS_UNF,XPRS_INF,XPRS_UNB,XPRS_OTH}) of string;!for status of given solution
    timetracker:    real; ! used to log timestamps for time consumption output
end-declarations

writeln("Building model...");
timetracker := timestamp; ! assigns current "timestamp" to timetracker

!Declarations of sets
declarations
    nMaxPeriods:        integer;
    nBidPoints:         integer;
end-declarations
    
initializations from Data_file
    nMaxPeriods;
    nBidPoints;
end-initializations
    
declarations
    Cycles:         set of integer;
    Periods:        set of integer;
    Turbines:       set of integer;
    FlowBreakPoints:set of integer;
    HeadBreakPoints:set of integer;
    DiagCons:       set of integer;
    Scenarios:      set of integer;
    BidPoints:      set of integer;!contains minDAPrice, 33%quartile, 67%quartile, maxDAprice+1
    CyclesReduced:  set of integer;
end-declarations    

Cycles          := 1 .. nCycles;!c
Periods         := 1 .. nMaxPeriods;!t
Turbines        := 1 .. nTurbines;!n
FlowBreakPoints := 1 .. nFlowBreakPoints;!r
HeadBreakPoints := 1 .. nHeadBreakPoints;!k
DiagCons        := 3 .. (nFlowBreakPoints+nHeadBreakPoints-1);!i    
Scenarios       := 1 .. nScenarios;!s
BidPoints       := 1 .. nBidPoints;!b

finalize(Cycles);
finalize(Periods);
finalize(Turbines);
finalize(FlowBreakPoints);
finalize(HeadBreakPoints);
finalize(DiagCons);
finalize(Scenarios);
finalize(BidPoints);

initializations from RedCycle_file
    CyclesReduced;
end-initializations

!Declarations of parameters in matrix form  
declarations
    nPeriods:   array(Cycles)           of real;    !number of production Periods in cycle c
    Flow:       array(FlowBreakPoints)  of real;!Flow in break point r
    PriceDA:    array(Cycles,Periods,Scenarios) of real;    !power price
    PriceRT:    array(Cycles,Periods,Scenarios) of real;    !power price
    BidPrice:   array(Cycles,Periods,BidPoints) of real;    !bid prices
    ResHead:    array(Cycles,Periods,HeadBreakPoints) of real;!Reservoir head in break point k
    Power:      array(Cycles, Periods, FlowBreakPoints, HeadBreakPoints )   of real;
    StartH:     array(Cycles)   of real;
    SOSindicesH:array(HeadBreakPoints) of integer;
    SOSindicesQ:array(FlowBreakPoints) of integer;
    MaxPower:   array(Cycles) of real;
end-declarations

initializations from Price_file
    PriceDA;
    PriceRT;
    BidPrice;
end-initializations
initializations from Cycle_file
    StartH;
    Flow;
    Power;
    ResHead;
    nPeriods;
end-initializations
initializations from Power_file
    MaxPower;
end-initializations

!Constants
declarations
    Duration:   integer;
    OprCost:    integer;
    Qmax:       integer;
    Qmin:       integer;
    Area:       real;
    Factor_Period: real;
    PowerFactor:integer;
    ResScale:   real;
    CycleScale: real;
    WeekScale:  integer;
end-declarations
    
initializations from Data_file
    Duration;
    OprCost;
    Qmax;
    Qmin;
    Area;
    Factor_Period;
    PowerFactor;
    ResScale;
    CycleScale;
    WeekScale;
    SOSindicesH;
    SOSindicesQ;
end-initializations

!Variables
declarations    
    spill:      dynamic array (CyclesReduced,Scenarios) of mpvar;
    bidDA:      dynamic array(CyclesReduced, Periods, BidPoints) of mpvar;
    soldDA:     dynamic array(CyclesReduced, Periods, Scenarios) of mpvar;
    soldRT:     dynamic array(CyclesReduced, Periods, Scenarios) of mpvar;
    power:      dynamic array (CyclesReduced,Periods,Turbines,Scenarios) of mpvar;
    resHeadTurbine: dynamic array (CyclesReduced,Periods,Turbines,Scenarios)    of mpvar;!turbine specific resHead
    resHead:    dynamic array (CyclesReduced,Periods,Scenarios) of mpvar;!turbine indep res head
    flow:       dynamic array (CyclesReduced,Periods,Turbines,Scenarios) of mpvar;
    starting:   dynamic array (CyclesReduced,Periods,Turbines,Scenarios) of mpvar;
    running:    dynamic array (CyclesReduced,Periods,Turbines,Scenarios) of mpvar;
    weight:     dynamic array (CyclesReduced,Periods,Turbines,Scenarios,FlowBreakPoints,HeadBreakPoints)    of mpvar;
    sumHWeight: dynamic array (CyclesReduced,Periods, Turbines,Scenarios,FlowBreakPoints)   of mpvar;
    sumQWeight: dynamic array (CyclesReduced,Periods, Turbines,Scenarios,HeadBreakPoints)   of mpvar;
    diagSum:    dynamic array (CyclesReduced,Periods, Turbines,Scenarios,DiagCons)  of mpvar;
end-declarations

forall (cc in CyclesReduced, ss in Scenarios) do
    create(spill(cc,ss));
end-do

forall (cc in CyclesReduced,tt in Periods|tt<=nPeriods(cc),nn in Turbines, ss in Scenarios, rr in FlowBreakPoints,kk in HeadBreakPoints|Power(cc,tt,rr,kk)>=0) do
    create(weight(cc,tt,nn,ss,rr,kk));  
end-do

forall(cc in CyclesReduced, tt in Periods | tt<=nPeriods(cc), bb in BidPoints) do 
    create(bidDA(cc,tt,bb));
end-do

forall (cc in CyclesReduced, tt in Periods | tt<=nPeriods(cc), ss in Scenarios) do
    create(soldDA(cc,tt,ss));
    create(soldRT(cc,tt,ss));
    create(resHead(cc,tt,ss));
    resHead(cc,tt,ss) is_free;
end-do

forall (cc in CyclesReduced, tt in Periods | tt<=nPeriods(cc), nn in Turbines, ss in Scenarios | (sum(rr in FlowBreakPoints, kk in HeadBreakPoints|exists(weight(cc,tt,nn,ss,rr,kk)))Power(cc,tt,rr,kk)>0)) do
    create(flow(cc,tt,nn,ss));
    create(resHeadTurbine(cc,tt,nn,ss));
    create(power(cc,tt,nn,ss));
    create(starting(cc,tt,nn,ss));
    create(running(cc,tt,nn,ss));
    resHeadTurbine(cc,tt,nn,ss) is_free;
    starting(cc,tt,nn,ss) is_binary;
    running(cc,tt,nn,ss) is_binary;
end-do

!variable holding sum of weighting variables over each grid column
forall(cc in CyclesReduced,tt in Periods|tt<=nPeriods(cc),nn in Turbines, ss in Scenarios,rr in FlowBreakPoints|sum(kk in HeadBreakPoints|exists(weight(cc,tt,nn,ss,rr,kk)))Power(cc,tt,rr,kk)>=0) do
    create(sumHWeight(cc,tt,nn,ss,rr));
end-do

!variable holding sum of weighting variables over each grid row
forall(cc in CyclesReduced,tt in Periods|tt<=nPeriods(cc),nn in Turbines, ss in Scenarios,kk in HeadBreakPoints|sum(rr in FlowBreakPoints|exists(weight(cc,tt,nn,ss,rr,kk)))Power(cc,tt,rr,kk)>=0) do
    create(sumQWeight(cc,tt,nn,ss,kk));
end-do

!variable holding sum of weighting variables over grid diagonals, used for ensuring unique interpolation
forall(cc in CyclesReduced,tt in Periods|tt<=nPeriods(cc),nn in Turbines, ss in Scenarios,ii in DiagCons) do
    create(diagSum(cc,tt,nn,ss,ii));
end-do

!Obj Function and restrictions
declarations
    ObjValue:       linctr;
    FlowCon:        dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
    ResHeadCon:     dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
    ResHeadFixCon:  dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
    PowerCon:       dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
    ProdSalesBal:   dynamic array (CyclesReduced, Periods, Scenarios) of linctr;
    SalesDA:        dynamic array (CyclesReduced, Periods, Scenarios) of linctr;
    SumWeight:      dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
    HeadFlowCon:    dynamic array (CyclesReduced, Periods, Scenarios) of linctr;
    HeadCon:        dynamic array (CyclesReduced, Periods, Scenarios) of linctr;
    TotalFlow:      dynamic array (CyclesReduced, Scenarios) of linctr;
    MaxFlow:        dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
    MinFlow:        dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
    TotalPower:     dynamic array (CyclesReduced, Scenarios) of linctr;
    StartCon:       dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
    Symmetry:       dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
    SumHeadWeight:  dynamic array (CyclesReduced, Periods, Turbines, Scenarios,FlowBreakPoints) of linctr;
    SumFlowWeight:  dynamic array (CyclesReduced, Periods, Turbines, Scenarios,HeadBreakPoints) of linctr;
    SumDiag:        dynamic array (CyclesReduced, Periods, Turbines, Scenarios,DiagCons) of linctr;
    SOS2R:          dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
    SOS2K:          dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
    SOS2Diag:       dynamic array (CyclesReduced, Periods, Turbines, Scenarios) of linctr;
end-declarations

!Annual Profit in GBP
ObjValue :=
    CycleScale*WeekScale*sum(ss in Scenarios)Prob*(
        sum(cc in CyclesReduced,tt in Periods)(PriceDA(cc,tt,ss)*soldDA(cc,tt,ss)+PriceRT(cc,tt,ss)*soldRT(cc,tt,ss))
        -
        sum(cc in CyclesReduced,tt in Periods,nn in Turbines) OprCost*starting(cc,tt,nn,ss)
        )/ResScale;

!For power opt
(!ObjValue :=
    sum(cc in CyclesReduced, tt in Periods, nn in Turbines, ss in Scenarios)Prob*power(cc,tt,nn,ss)/ResScale;!)

!Flow is a weighted sum of flow grid point values
forall(cc in CyclesReduced,tt in Periods| tt<=nPeriods(cc),nn in Turbines,ss in Scenarios | exists(flow(cc,tt,nn,ss))) do
    FlowCon(cc,tt,nn,ss):=
        flow(cc,tt,nn,ss)=sum(rr in FlowBreakPoints,kk in HeadBreakPoints) Flow(rr)*weight(cc,tt,nn,ss,rr,kk);
end-do

!Reservoir head is a weighted sum of Reservoir head grid point values
forall(cc in CyclesReduced, tt in Periods| tt<=nPeriods(cc), nn in Turbines,ss in Scenarios | exists(resHeadTurbine(cc,tt,nn,ss))) do
    ResHeadCon(cc,tt,nn,ss):=
        resHeadTurbine(cc,tt,nn,ss)=sum(rr in FlowBreakPoints,kk in HeadBreakPoints) ResHead(cc,tt,kk)*weight(cc,tt,nn,ss,rr,kk);
end-do

!Turbine dependent reservoir head variables must be equal
forall(cc in CyclesReduced, tt in Periods| tt<=nPeriods(cc), nn in Turbines,ss in Scenarios | exists(resHeadTurbine(cc,tt,nn,ss))) do
    ResHeadFixCon(cc,tt,nn,ss):=
        resHead(cc,tt,ss)=resHeadTurbine(cc,tt,nn,ss);
end-do

!power generated is a weighted sum of power grid point values, in MW
forall(cc in CyclesReduced, tt in Periods| tt<=nPeriods(cc), nn in Turbines,ss in Scenarios | exists(power(cc,tt,nn,ss))) do
    PowerCon(cc,tt,nn,ss):=
        power(cc,tt,nn,ss)=sum(rr in FlowBreakPoints,kk in HeadBreakPoints) 0.9*Power(cc,tt,rr,kk)*weight(cc,tt,nn,ss,rr,kk)/PowerFactor;
end-do

!upper limit on total power generated over a cycle
forall(cc in CyclesReduced, ss in Scenarios) do
    TotalPower(cc,ss):=
        sum(tt in Periods, nn in Turbines) power(cc,tt,nn,ss)<=MaxPower(cc)*ResScale;
end-do

!Power generated is scaled to MW and sold in the day-ahead and real-time market
forall(cc in CyclesReduced, tt in Periods| tt<=nPeriods(cc),ss in Scenarios) do
    ProdSalesBal(cc,tt,ss):=
        soldDA(cc,tt,ss) + soldRT(cc,tt,ss) = sum(nn in Turbines)power(cc,tt,nn,ss)*Factor_Period;
end-do

!Delivery obligation in day-ahead market based on bidding volumes and market realisations
!Obligation calculated as interpolation between bidding volumes
forall(cc in CyclesReduced, tt in Periods| tt<=nPeriods(cc),ss in Scenarios, bb in BidPoints|bb>1) do
    if (PriceDA(cc,tt,ss)>=BidPrice(cc,tt,bb-1) and PriceDA(cc,tt,ss)<BidPrice(cc,tt,bb)) then
        SalesDA(cc,tt,ss):=
            soldDA(cc,tt,ss)=(PriceDA(cc,tt,ss)-BidPrice(cc,tt,bb-1))/(BidPrice(cc,tt,bb)-BidPrice(cc,tt,bb-1))*bidDA(cc,tt,bb)
            +
            (BidPrice(cc,tt,bb)-PriceDA(cc,tt,ss))/(BidPrice(cc,tt,bb)-BidPrice(cc,tt,bb-1))*bidDA(cc,tt,bb-1);
    end-if
end-do!)

!Weighting variables sum to 1, grid includes zero flow
forall(cc in CyclesReduced,tt in Periods| tt<=nPeriods(cc), nn in Turbines,ss in Scenarios) do
    SumWeight(cc,tt,nn,ss):=
        sum(rr in FlowBreakPoints,kk in HeadBreakPoints) weight(cc,tt,nn,ss,rr,kk)=1;!Zero flow included in Break Points
end-do

!reservoir head decreases/increases with flow out of/into reservoir
forall(cc in CyclesReduced, tt in Periods|tt>1 and tt<=nPeriods(cc),ss in Scenarios) do
    if StartH(cc)<0 then 
        HeadFlowCon(cc,tt,ss):=
            resHead(cc,tt,ss)=StartH(cc)+Duration/(Area*ResScale)*sum(ii in 1..(tt-1),nn in Turbines)flow(cc,ii,nn,ss);
    else 
        HeadFlowCon(cc,tt,ss):=
            resHead(cc,tt,ss)=StartH(cc)-Duration/(Area*ResScale)*sum(ii in 1..(tt-1),nn in Turbines)flow(cc,ii,nn,ss);     
    end-if
end-do!)

!Reservoir head equals tide level for t=1
forall(cc in CyclesReduced, ss in Scenarios) do
    HeadFlowCon(cc,1,ss):=
        resHead(cc,1,ss)=StartH(cc);
end-do      

!Limiting total flow in a cycle to reservoir volume = HeadDifference*ReservoirArea. 
!ResScale reduces total reservoir volume to facilitate reduced number of turbines
forall(cc in CyclesReduced,ss in Scenarios) do
    if cc<nCycles then
        if StartH(cc)<0 then
            TotalFlow(cc,ss):=
                sum(tt in Periods,nn in Turbines)flow(cc,tt,nn,ss)*Duration + spill(cc,ss)= -(StartH(cc)-StartH(cc+1))*Area*ResScale-70200;
        else
            TotalFlow(cc,ss):=
                sum(tt in Periods,nn in Turbines)flow(cc,tt,nn,ss)*Duration + spill(cc,ss)= +(StartH(cc)-StartH(cc+1))*Area*ResScale-70200;
        end-if
    else
        if StartH(cc)<0 then
            TotalFlow(cc,ss):=
                sum(tt in Periods,nn in Turbines)flow(cc,tt,nn,ss)*Duration + spill(cc,ss)= -2*StartH(cc)*Area*ResScale-70200;
        else
            TotalFlow(cc,ss):=
                sum(tt in Periods,nn in Turbines)flow(cc,tt,nn,ss)*Duration + spill(cc,ss)= +2*StartH(cc)*Area*ResScale-70200;
        end-if          
    end-if              
end-do

!max flow
forall(cc in CyclesReduced,tt in Periods | tt<=nPeriods(cc),nn in Turbines,ss in Scenarios | exists(flow(cc,tt,nn,ss))) do
    MaxFlow(cc,tt,nn,ss):=
        flow(cc,tt,nn,ss)-Qmax*running(cc,tt,nn,ss) <= 0;
end-do

!Min flow
forall(cc in CyclesReduced,tt in Periods | tt<=nPeriods(cc),nn in Turbines,ss in Scenarios | exists(flow(cc,tt,nn,ss))) do
    MinFlow(cc,tt,nn,ss):=
        flow(cc,tt,nn,ss)-Qmin*running(cc,tt,nn,ss) >= 0;
end-do

!control of binary starting and running variables in time periods > 1
forall(cc in CyclesReduced,tt in Periods | tt>1 and tt<=nPeriods(cc),nn in Turbines,ss in Scenarios | exists(running(cc,tt,nn,ss))) do
    StartCon(cc,tt,nn,ss):=
        running(cc,(tt-1),nn,ss)+starting(cc,tt,nn,ss)-running(cc,tt,nn,ss)>=0;
end-do

!control of binary starting and running variables in the first time period
forall(cc in CyclesReduced,nn in Turbines,ss in Scenarios | exists(running(cc,1,nn,ss))) do
    StartCon(cc,1,nn,ss):=
        starting(cc,1,nn,ss)=running(cc,1,nn,ss);
end-do

!symmetry breaking constraint for turbines
forall(cc in CyclesReduced,tt in Periods | tt<=nPeriods(cc),nn in Turbines|nn<nTurbines,ss in Scenarios | exists(flow(cc,tt,nn,ss)) and exists(flow(cc,tt,nn+1,ss))) do
    Symmetry(cc,tt,nn,ss):=
        flow(cc,tt,nn,ss)-flow(cc,tt,(nn+1),ss)>=0;
end-do

!sum of weighting variables over each grid column
forall(cc in CyclesReduced,tt in Periods| tt<=nPeriods(cc),nn in Turbines,ss in Scenarios,rr in FlowBreakPoints | exists(sumHWeight(cc,tt,nn,ss,rr))) do
    SumHeadWeight(cc,tt,nn,ss,rr):=
        sumHWeight(cc,tt,nn,ss,rr) = sum(kk in HeadBreakPoints) weight(cc,tt,nn,ss,rr,kk);
end-do

!sum of weighting variables over each grid row
forall(cc in CyclesReduced,tt in Periods| tt<=nPeriods(cc),nn in Turbines,ss in Scenarios, kk in HeadBreakPoints | exists(sumQWeight(cc,tt,nn,ss,kk))) do
    SumFlowWeight(cc,tt,nn,ss,kk):=
        sumQWeight(cc,tt,nn,ss,kk) = sum(rr in FlowBreakPoints) weight(cc,tt,nn,ss,rr,kk);
end-do

!sum of weighting variables over each grid diagonal
forall(cc in CyclesReduced,tt in Periods| tt<=nPeriods(cc),nn in Turbines,ss in Scenarios, ii in DiagCons) do
    SumDiag(cc,tt,nn,ss,ii):=
        diagSum(cc,tt,nn,ss,ii) = sum(rr in FlowBreakPoints|exists(weight(cc,tt,nn,ss,rr,(ii+rr-nFlowBreakPoints))))weight(cc,tt,nn,ss,rr,(ii+rr-nFlowBreakPoints));
end-do!)

!SOS2 imposed on sum of weighting variables over each grid column
forall(cc in CyclesReduced, tt in Periods| tt<=nPeriods(cc), nn in Turbines,ss in Scenarios) do
    SOS2R(cc,tt,nn,ss):=
        sum(rr in FlowBreakPoints) SOSindicesQ(rr)*sumHWeight(cc,tt,nn,ss,rr) is_sos2;
end-do

!SOS2 imposed on sum of weighting variables over each grid row
forall(cc in CyclesReduced, tt in Periods| tt<=nPeriods(cc), nn in Turbines,ss in Scenarios) do
    SOS2K(cc,tt,nn,ss):=
        sum(kk in HeadBreakPoints) SOSindicesH(kk)*sumQWeight(cc,tt,nn,ss,kk) is_sos2;
end-do

!SOS2 imposed on sum of weighting variables over each grid diagonal
forall(cc in CyclesReduced, tt in Periods| tt<=nPeriods(cc), nn in Turbines,ss in Scenarios) do
    SOS2Diag(cc,tt,nn,ss):=
        sum(ii in DiagCons)ii*diagSum(cc,tt,nn,ss,ii) is_sos2;
end-do!)

writeln("\nModel building completed in ", timestamp - timetracker, " seconds");
writeln("\nSolving model...");
timetracker := timestamp;
        

status::([XPRS_OPT,XPRS_UNF,XPRS_INF,XPRS_UNB,XPRS_OTH])[
                "Optimum found","Unfinished","Infeasible","Unbounded","Failed"];

maximize(ObjValue);

writeln("\nModel solved in ", timestamp - timetracker," seconds");
writeln("Status: " + status(getprobstat));

!write problem overview to file
fopen(ProblemFile, F_APPEND);
if parseint(nameapp,1) = 1 then
    writeln("Price Case: ", PriceCase, " Markets: ", Markets);
    writeln("Scenarios: ", nScenarios);
    writeln("Turbines: ", nTurbines);
    writeln("Time Limit: ", MaxTime);
end-if
writeln("---------------------------------------------------------------------------");
forall (cc in CyclesReduced) do
    writeln("Cycle: ", cc);
end-do
writeln("Profit = ", getobjval);
writeln("Status: ", + status(getprobstat));
writeln("Model solved in ", timestamp - timetracker," seconds");
if getobjval<=0 then
    writeln("Zero-solution or no solution found, gap is not available");
else
    writeln("Gap: ", (getparam("xprs_bestbound")-getobjval)/getobjval*100, " %");!gap in %  
end-if
writeln("Upper bound: ", getparam("xprs_bestbound"));
fclose(F_APPEND);

!write profits to file
fopen(ProfitFile, F_APPEND);
if parseint(nameapp,1) = 1 then
    writeln(strfmt("Cycle",6), strfmt("Profits", 12));
end-if
forall(cc in CyclesReduced) do
    writeln(strfmt(cc, 6), strfmt(getobjval, 12,2));
end-do
fclose(F_APPEND);

!write generation schedule to file
fopen(PowerFile, F_APPEND);
if parseint(nameapp,1) = 1 then
    writeln(strfmt("Cycle",6), strfmt("Period",7), strfmt("Scenario",10), strfmt("Power Gen",12), strfmt("#Turbines",12));
end-if
forall(cc in CyclesReduced, tt in Periods|tt<=nPeriods(cc), ss in Scenarios) do
    writeln(strfmt(cc,6), strfmt(tt,7), strfmt(ss,10), strfmt(sum(nn in Turbines)getsol(power(cc,tt,nn,ss))/ResScale,12,2), 
    strfmt(sum(nn in Turbines)getsol(running(cc,tt,nn,ss)),12));
end-do
fclose(F_APPEND);

!write sales solutions to file
fopen(SalesFile, F_APPEND);
if parseint(nameapp,1) = 1 then
    writeln(strfmt("Cycle",6), strfmt("Period", 7), strfmt("Scenario", 10), strfmt("SoldDA",12),
    strfmt("SoldRT",12));
end-if
forall(cc in CyclesReduced, tt in Periods|tt<=nPeriods(cc), ss in Scenarios) do
    writeln(strfmt(cc,6), strfmt(tt,7), strfmt(ss,10), strfmt(getsol(soldDA(cc,tt,ss))/ResScale,12,2), 
    strfmt(getsol(soldRT(cc,tt,ss))/ResScale,12,2));
end-do
fclose(F_APPEND);

!write bids to file
fopen(BidFile, F_APPEND);
if parseint(nameapp,1) = 1 then
    writeln(strfmt("Cycle",6), strfmt("Period", 7), strfmt("BidPrice", 10), strfmt("Bid Volumes",12));
end-if
forall(cc in CyclesReduced, tt in Periods|tt<=nPeriods(cc), bb in BidPoints) do
    writeln(strfmt(cc,6), strfmt(tt,7), strfmt(bb,10), strfmt(getsol(bidDA(cc,tt,bb))/ResScale,12,2));
end-do
fclose(F_APPEND);

end-model