(!****************************************************** Mosel Example Problems ====================== file a6electr_ro_hdem.mos `````````````````````````` Security constrained robust unit commitment under load variation. The present model produces a unit commitment schedule that does not require startup or shutdown change as long as the realized load is a convex combination of some historical load expressed with the 'scenario' uncertainty set. Features: - Using the 'scenario' uncertainty set - Unit commitment robust against load variaton Objectives: - Minimize expected running cost - No commitment modification if load varies What to look at: - The number of committed units per periods - The power generation reserve - The loss-of-load risk source: see a6electr.mos for the original formulation (c) 2014 Fair Isaac Corporation author: S. Heipcke, Mar. 2002 S. Lannez, Mar. 2014 *******************************************************!) model "A-6 Electricity production (robust)" uses "mmsystem", "mmrobust" declarations NT = 7 ! Number of time periods TIME = 1..NT ! Time periods PLANTS = 1..4 ! Power generation plant UNITS: set of string ! Power generation units LEN, DEM: array(TIME) of integer ! Length and demand of time periods PMIN,PMAX: array(PLANTS) of integer ! Min. and max output of a generator type CSTART: array(PLANTS) of integer ! Start-up cost of a generator CMIN: array(PLANTS) of integer ! Hourly cost of gen. at min. output CADD: array(PLANTS) of real ! Cost/hour/MW of prod. above min. level PLANT: array(UNITS) of integer ! Unit generation plant YEARS: range ! Historical years HDEM: array(YEARS,TIME) of integer ! Historical demand profiles start: array(UNITS,TIME) of mpvar ! Is generation unit starting ? work: array(UNITS,TIME) of mpvar ! Is generation unit up ? padd: array(UNITS,TIME) of mpvar ! Production above min. output level end-declarations initializations from 'a6electr_ro_simple.dat' LEN DEM PMIN PMAX CSTART CMIN CADD PLANT HDEM end-initializations write("\n", strfmt("Unit type",-20)) forall(p in PLANTS) write(strfmt(p,8)," ") write("\n", strfmt("Pmax",20)) forall(p in PLANTS) write(strfmt(PMAX(p),8)," ") writeln writeln !***************************Subroutines*************************** !**** Helper routine to create scenario uncertain set **** public procedure makescenario(a:array(r1:range,c1:range) of integer, u:array(c2:range) of uncertain) declarations aa:dynamic array(r0:range,U:set of uncertain) of real end-declarations forall(i in r1, j in c1|exists(a(i,j)) and exists(u(j)),n as counter) do aa(n,u(j)):=a(i,j) end-do scenario(aa) end-procedure !**** Robust against past scenario **** public procedure robust_demand declarations demand: array(TIME) of uncertain ! Uncertain power demand end-declarations ! Add uncertainty set (time correlation) makescenario(HDEM,demand) ! Security reserve upward forall(t in TIME) sum(u in UNITS) PMAX(PLANT(u))*work(u,t) >= demand(t) end-procedure !**** Print highest loss of load in case of worst case scenario realization procedure print_risk write("\n", strfmt("Loss of load",20)) forall(t in TIME) do s:= sum(u in UNITS) (work(u,t).sol*PMAX(PLANT(u))) ! Total capacity v:= maxlist(DEM(t), max(y in YEARS) HDEM(y,t)) ! Hist. peak demand if(s < v) then write(strfmt(v-s,8)) else write(strfmt("",8)) end-if end-do end-procedure !**** Solution printing **** procedure print_solution writeln("Total cost: ", getobjval) write(strfmt("Time period ",-20)) ct:=0 forall(t in TIME) do write(strfmt(ct,5), "-", strfmt(ct+LEN(t),2), "") ct+=LEN(t) end-do write("\n",strfmt("Up Units",-20)) forall(t in TIME) write(strfmt(sum(u in UNITS) work(u,t).sol,8)) write("\n", strfmt("Gen. Pwr.",20)) forall(t in TIME) write(strfmt(sum(u in UNITS) (work(u,t).sol*PMIN(PLANT(u))+padd(u,t).sol),8)) write("\n", strfmt("Gen. Cap.",20)) forall(t in TIME) write(strfmt(sum(u in UNITS) (work(u,t).sol*PMAX(PLANT(u))),8)) write("\n", strfmt("Res. Dn",20)) forall(t in TIME) write(strfmt(sum(u in UNITS) (padd(u,t).sol),8)) write("\n", strfmt("Res. Up",20)) forall(t in TIME) write(strfmt(sum(u in UNITS) work(u,t).sol*PMAX(PLANT(u)) - DEM(t),8)) print_risk writeln end-procedure !***************************Main model*************************** ! Create decision variables forall(u in UNITS, t in TIME) do start(u,t) is_binary work(u,t) is_binary end-do ! Objective function: total daily cost Cost:= sum(u in UNITS, t in TIME) (CSTART(PLANT(u))*start(u,t) + LEN(t)*(CMIN(PLANT(u))*work(u,t) + CADD(PLANT(u))*padd(u,t))) ! Number of generators started per period and per type forall(u in UNITS, t in TIME) start(u,t) >= work(u,t) - if(t>1, work(u,t-1), work(u,NT)) ! Limit on power production range forall(u in UNITS, t in TIME) padd(u,t) <= (PMAX(PLANT(u))-PMIN(PLANT(u)))*work(u,t) ! Power balance forall(t in TIME) sum(u in UNITS) (PMIN(PLANT(u))*work(u,t) + padd(u,t)) = DEM(t) ! Symmetry breaker forall(t in TIME, p in PLANTS) forall(u1, u2 in UNITS | PLANT(u1) = PLANT(u2) and u1= work(u2,t) !**** Solving and solution reporting **** ! Solve the original problem minimize(Cost) writeln("\n=== Nominal case ===\n") print_solution ! Add robust constraints robust_demand ! Solve robust the problem setparam("ROBUST_UNCERTAIN_OVERLAP", true) minimize(Cost) writeln("\n=== Robust against demand variation ===\n") print_solution end-model