(!******************************************************
   Mosel Example Problems
   ======================

   file a6electrg.mos
   ``````````````````
   Production of electricity.
   Plotting parametrics in the reserve capacity.
   
   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, Oct. 2004, rev. Sep. 2017
*******************************************************!)

model "A-6 Electricity production (with graph drawing)"
 uses "mmxprs", "mmsvg"
 
 declarations
  NT = 7
  TIME = 1..NT                       ! Time periods
  TYPES = 1..4                       ! Power generator types

  LEN, DEM: array(TIME) of integer   ! Length and demand of time periods
  PMIN,PMAX: array(TYPES) of integer ! Min. & max output of a generator type
  CSTART: array(TYPES) of integer    ! Start-up cost of a generator
  CMIN: array(TYPES) of integer      ! Hourly cost of gen. at min. output
  CADD: array(TYPES) of real         ! Cost/hour/MW of prod. above min. level
  AVAIL: array(TYPES) of integer     ! Number of generators per type
  bas: basis                         ! LP basis

  start: array(TYPES,TIME) of mpvar  ! No. of gen.s started in a period
  work: array(TYPES,TIME) of mpvar   ! No. of gen.s working during a period
  padd: array(TYPES,TIME) of mpvar   ! Production above min. output level
 end-declarations
 
! Time period data
 LEN:: [    6,     3,     3,     2,     4,     4,     2]
 DEM:: [12000, 32000, 25000, 36000, 25000, 30000, 18000]

! Power plant data
 PMIN::   [ 750, 1000, 1200, 1800]
 PMAX::   [1750, 1500, 2000, 3500]
 CSTART:: [5000, 1600, 2400, 1200]
 CMIN::   [2250, 1800, 3750, 4800]
 CADD::   [ 2.7,  2.2,  1.8,  3.8]
 AVAIL::  [  10,    4,    8,    3]
 
 
! Objective function: total daily cost
 Cost:= sum(p in TYPES, t in TIME) (CSTART(p)*start(p,t) +
          LEN(t)*(CMIN(p)*work(p,t) + CADD(p)*padd(p,t))) 
                                   
! Number of generators started per period and per type
 forall(p in TYPES, t in TIME) 
  NumStart(p,t):= start(p,t) >= work(p,t) - if(t>1, work(p,t-1), work(p,NT))

! Limit on power production above minimum level
 forall(p in TYPES, t in TIME) 
  Limit(p,t):= padd(p,t) <= (PMAX(p)-PMIN(p))*work(p,t)

! Satisfy demands
 forall(t in TIME) 
  Demand(t):= sum(p in TYPES) (PMIN(p)*work(p,t) + padd(p,t)) >= DEM(t)

! Security reserve (reserve initially at 0%)
 forall(t in TIME) 
  Reserve(t):= sum(p in TYPES) PMAX(p)*work(p,t) >= DEM(t)

! Limit number of available generators; numbers of generators are integer
 forall(p in TYPES, t in TIME) do
  work(p,t) <= AVAIL(p)
  work(p,t) is_integer
 end-do


! Uncomment the following line to see the Optimizer log
! setparam("XPRS_verbose", true)

! Solve as an LP problem and save the basis
 setparam("XPRS_PRESOLVE", 0)
 minimize(XPRS_LPSTOP, Cost)
 savebasis(bas)
 setparam("XPRS_PRESOLVE", 1)

! Define a user graph
 svgaddgroup("Graph", "Reserve %")

! Solve the problem with different values of reserve (from 0% to 20%)
 forall(r in 0..20) do
  forall(t in TIME)                   ! Redefine the reserve constraints
   Reserve(t):= sum(p in TYPES) PMAX(p)*work(p,t) >= (1+r/100)*DEM(t)

  loadprob(Cost)                      ! Load problem into the Optimizer
  loadbasis(bas)                      ! Load the saved basis
  minimize(Cost)                      ! Solve the modified problem

  writeln(r, "%: ", getobjval)        ! Print solution value
  svgaddpoint("Graph", r, getobjval/1000)  ! Display solution in user graph
 end-do

! Scale the size of the displayed graph
 svgsetgraphscale(10)
 svgsetgraphpointsize(3)
 svgsetgraphlabels("Reserve %", "Cost")

 svgsave("a6electr.svg")
 svgrefresh
 svgwaitclose("Close browser window to terminate model execution.", 1)
end-model