Implementation
The Original Unit Commitment Implementation
The Mosel code printed below implements and solves the original formulation of the unit commitment problem.
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
! 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<u2)
work(u1,t) >= 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
The following output printing routines are invoked from the model shown above:
!**** 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
The Load Robust Unit Commitment Implementation
The extension to scenario-based robust optimization is implemented by the following 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
The N-k Contingency-Constrained Unit Commitment Implementation
Replacing the call to the procedure robust_demand by a call to the procedure robust_contingency printed below will turn the problem formulation into a N-k contingency-constrained unit commitment model.
The example also demonstrates how to use the value of the uncertain to extract meaningful information about the units which are the most critical ones for the system.
!**** Ensure enough power production capacity allows to
!**** satisfy load even if 'k' units are forced in outage.
procedure robust_contingency(k: integer)
forall(t in TIME, u in UNITS) create(outage(u,t))
forall(t in TIME, u in UNITS) outage(u,t) >= 0
! If the unit is forced in outage, the total capacity
! of the generating units are lost
forall(t in TIME, u in UNITS) outage(u,t) <= 1
! forall(t in TIME) sum(u in UNITS) outage(u,t) <= k
forall(t in TIME) cardinality(union(u in UNITS | exists(outage(u,t))) {outage(u,t)}, k)
forall(t in TIME)
RobCtrMaxDem(t) :=
sum(u in UNITS) PMAX(PLANT(u))*work(u,t) -
sum(u in UNITS) PMAX(PLANT(u))*work(u,t)*outage(u,t) >= DEM(t)
! In order to reduce the number of variables we reuse
! outage() and let the solver do the variable duplication
! by setting the ROBUST_UNCERTAIN_OVERLAP parameter
setparam("ROBUST_UNCERTAIN_OVERLAP",true)
forall(t in TIME)
RobCtrMinDem(t) :=
sum(u in UNITS) PMIN(PLANT(u))*work(u,t) -
sum(u in UNITS) PMIN(PLANT(u))*work(u,t)*outage(u,t) <= DEM(t)
end-procedure