######################################################################## # This code is part of UCP_VALVE tube packing software # # Copyright (C) 2013 # JP Pedroso(*), M Kubo, A Viana # (*) Universidade do Porto, Faculdade de Ciencias # E-mail: . # # UCP_VALVE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published # by the Free Software Foundation, either version 3 of the License, # or (at your option) any later version. # # UCP_VALVE is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with UCP_VALVE. If not, see . ######################################################################## """ ucp_valve.py: python code for thermal ucp/load dispatch with valve effect Based on: [1] J.~P. Pedroso, A.~Viana, and M.~Kubo. Unit commitment with valve-point loading effect. (Submitted to publication). USAGE: python ucp_valve.py CPU inst nperiods where: - CPU is the limit CPU allowed (s) - inst is one of: toy eld13 eld40 ucp10 ucp13 ucp40 - nperiods is the number of periods to consider, from 1 to 24 REQUIRES: a version of gurobipy installed [see http://www.gurobi.com] """ import math from ucp_data import * from gurobipy import * def cost(a,b,c,d,e,f,p_min,p): """cost(a,b,c,d,e,f,p_min,p): operating cost taking into accept valve-point effect""" return a + b*p + c*p*p + d*p**3 + abs(e*math.sin(f*(p_min-p))) def dcost_dp(a,b,c,d,e,f,p_min,p): """dcost_dp: cost function's derivative""" if math.sin(f*(p_min-p)) > 0: return b + 2*c*p + 3*d*p**2 - (f*e*math.cos(f*(p_min-p))) else: return b + 2*c*p + 3*d*p**2 + (f*e*math.cos(f*(p_min-p))) def lower_brkpts(a,b,c,d,e,f,p_min,p_max,n): """lower_brkpts: lower approximation of the cost function""" EPS = 1.e-12 # for assuring left/right derivatives are finite if f == 0: f = math.pi/(p_max-p_min) brk = [] nvalve = int(math.ceil(f*(p_max-p_min)/math.pi)) for i in range(nvalve+1): p0 = p_min + i*math.pi/f if p0 >= p_max-EPS: brk.append((p_max, cost(a,b,c,d,e,f,p_min,p_max))) break for j in range(n): p = p0 + j*math.pi/f/n if p >= p_max: break brk.append((p, cost(a,b,c,d,e,f,p_min,p))) return brk def has_added_p(p_min,p_max,f,prev,p): """has_added_p(p_min,p_max,f,prev,p): check if p is alread "covered" by some value in prev""" # print "has_added_p: checking", p , "in", prev nvalve = int(math.ceil(f*(p_max-p_min)/math.pi)) for i in range(nvalve+1): left = p_min + i*math.pi/f right = p_min + (i+1)*math.pi/f if p > left+EPS and p < right-EPS: for r in prev: if r > p_min + i*math.pi/f + EPS and r < p_min + (i+1)*math.pi/f - EPS: # print "point covered by", r return False elif p >= left-EPS and p <= right+EPS: # valve point # print "valve point" return False return True def expand_brkpts(a,b,c,d,e,f,p_min,p_max,p_list,n): """expand_brkpts(a,b,c,d,e,f,p_min,p_max,p_list,n): lower approximation of the cost function""" EPS = 1.e-12 # for assuring left/right derivatives if f == 0: f = math.pi/(p_max-p_min) brk = [] nvalve = int(math.ceil(f*(p_max-p_min)/math.pi)) k = 0 if k < len(p_list): p_sol = p_list[k] else: p_sol = float("inf") for i in range(nvalve+1): p0 = p_min + i*math.pi/f if p0 >= p_max-EPS: brk.append((p_max, cost(a,b,c,d,e,f,p_min,p_max))) break for j in range(n): p = p0 + j*math.pi/f/n if p >= p_max-EPS: break brk.append((p, cost(a,b,c,d,e,f,p_min,p))) if p_min + i*math.pi/f >= p_sol-EPS or p_min + (i+1)*math.pi/f <= p_sol+EPS: break while p_sol < p_min + (i+1)*math.pi/f + EPS: k += 1 if k < len(p_list): p_sol = p_list[k] else: p_sol = float("inf") return brk def upper_brkpts(a,b,c,d,e,f,p_min,p_max,n): """upper_brkpts(a,b,c,d,e,f,p_min,p_max,n): upper approximation of the cost function""" EPS = 1.e-6 # for assuring left/right derivatives if f == 0: f = math.pi/(p_max-p_min) brk = [] nvalve = int(math.ceil(f*(p_max-p_min)/math.pi)) for i in range(nvalve+1): p0 = p_min + i*math.pi/f if p0 >= p_max-EPS: brk.append((p_max, cost(a,b,c,d,e,f,p_min,p_max))) break p2 = p0 brk.append((p2, cost(a,b,c,d,e,f,p_min,p2))) for j in range(1,n+1): p1 = p2 p2 = min(p_max,p0 + j*math.pi/f/n) if p1 >= p_max-EPS: break cost1,cost2 = cost(a,b,c,d,e,f,p_min,p1), cost(a,b,c,d,e,f,p_min,p2) dcost1,dcost2 = dcost_dp(a,b,c,d,e,f,p_min,p1+EPS), dcost_dp(a,b,c,d,e,f,p_min,p2-EPS) newp = (cost2 - cost1 + dcost1*p1 - dcost2*p2)/(dcost1 - dcost2) newcost = cost(a,b,c,d,e,f,p_min,p1) + dcost_dp(a,b,c,d,e,f,p_min,p1+EPS)*(newp-p1) brk.append((newp, newcost)) return brk def ucp_base(U, T, p_min, p_max, y_prev, t_cold, n_init, min_on, min_off, a_hot, a_cold, dem, res): """ucp_base(...): prepare mip with standard ucp model (kazarlis), not setting objective""" y, on, off, cold, hot, p, fuel, startup = {}, {}, {}, {}, {}, {}, {}, {} model = Model("Unit Commitment Problem") for u in U: y[u,0] = model.addVar(vtype="B", name="y(%s,%s)"%(u,0)) for t in range(1,T+1): y[u,t] = model.addVar(vtype="B", name="y(%s,%s)"%(u,t)) on[u,t] = model.addVar(lb=0, ub=1, name="on(%s,%s)"%(u,t)) off[u,t] = model.addVar(lb=0, ub=1, name="off(%s,%s)"%(u,t)) cold[u,t] = model.addVar(lb=0, ub=1, name="cold(%s,%s)"%(u,t)) hot[u,t] = model.addVar(lb=0, ub=1, name="hot(%s,%s)"%(u,t)) p[u,t] = model.addVar(lb=0,name="p(%s,%s)"%(u,t)) fuel[u,t] = model.addVar(lb=0,name="fuel(%s,%s)"%(u,t)) startup[u,t] = model.addVar(lb=0,name="startup(%s,%s)"%(u,t)) model.update() for t in range(1,T+1): # demand satisfaction model.addConstr(quicksum(p[u,t] for u in U) == dem[t], "demand(%s)"%t) model.addConstr(quicksum(y[u,t]*p_max[u] for u in U) >= dem[t]+res[t], "reserve(%s)"%t) # capacity for u in U: model.addConstr(y[u,t] * p_min[u] <= p[u,t], "min_prod(%s,%s)"%(u,t)) model.addConstr(y[u,t] * p_max[u] >= p[u,t], "max_prod(%s,%s)"%(u,t)) for u in U: if y_prev[u] == 1: cold[u,1].ub = 0 hot[u,1].ub = 0; for t in range(max(0,min_on[u]-n_init[u])+1): y[u,t].lb = 1 # for all units with t_cold > min_down and y_prev[u] = 1, # start up costs are never cold for t = 1...t_cold+1 if t_cold[u] > min_off[u]: for t in range(1,min(T, t_cold[u]+1) + 1): cold[u,t].ub = 0 else: # y_prev[u] == 0: if n_init[u] > t_cold[u]: model.addConstr(cold[u,1] == y[u,1]) hot[u,1].ub = 0 else: cold[u,1].ub = 0 model.addConstr(hot[u,1] == y[u,1]) for t in range(max(0,min_off[u]-n_init[u])+1): y[u,t].ub = 0; # for all units with t_cold > min_down time and y_prev[u] = 0 # start ups can be cold and hot for t = 1... t_cold+1 if t_cold[u] > min_off[u]: for t in range(2,min(T,t_cold[u]+1)+1): n = max(0, t_cold[u] - (t-1) - n_init[u] + 1) model.addConstr(y[u,t] - quicksum(y[u,i] for i in range(1,t)) - n*(1-y_prev[u]) <= cold[u,t], "cold_init(%s,%s)" % (u,t)) if t_cold[u] > 0: for t in range(t_cold[u]+2,T+1): model.addConstr(y[u,t] - quicksum(y[u,i] for i in range(t-t_cold[u]-1, t)) <= cold[u,t], "c_cold(%s,%s)"%(u,t)) for t in range(1,T+1): # for all units with t_cold < min_down time start ups are always cold if t_cold[u] < min_off[u] or t_cold[u] == 0: hot[u,t].ub = 0 # compute aux vars model.addConstr(y[u,t] - y[u,t-1] <= on[u,t], "c_on(%s,%s)" % (u,t)) model.addConstr(off[u,t] == on[u,t] + y[u,t-1] - y[u,t], "c_off(%s,%s)" % (u,t)) t0 = max(t-min_on[u]+1, 1) model.addConstr(quicksum(on[u,i] for i in range(t0,t+1)) <= y[u,t], "keep_on(%s,%s)" % (u,t)) t0 = max(t-min_off[u]+1, 1) model.addConstr(quicksum(off[u,i] for i in range(t0,t+1)) <= 1-y[u,t], "keep_off(%s,%s)" % (u,t)) model.addConstr(hot[u,t] + cold[u,t] == on[u,t], "hot_or_cold(%s,%s)" % (u,t)) model.__data = y, p, hot, cold, on, off, fuel, startup return model def ucp_quad_obj(model, U, T, a, b, c, a_hot, a_cold, p_min, p_max): """ucp_quad_obj(model, ...): set quadratic objective to given base model""" # set objective y, p, hot, cold, on, off, fuel, startup = model.__data obj = 0 for u in U: for t in range(1,T+1): model.addConstr(fuel[u,t] >= a[u]*y[u,t] + b[u]*p[u,t] + c[u]*p[u,t]*p[u,t]) # comment for quad.cost model.addConstr(startup[u,t] == a_hot[u]*hot[u,t] + a_cold[u]*cold[u,t]) obj += fuel[u,t] + startup[u,t] model.setObjective(obj, GRB.MINIMIZE) def ucp_lin_obj(model, U, T, brk): """ucp_lin_obj(model, ...): set objective based on piecewise linear approximation to given base model""" y, p, hot, cold, on, off, fuel, startup = model.__data obj = 0 for u in U: for t in range(1,T+1): # for k,(X,Y) in enumerate(brk[u,t]): # print k,X,Y abrk = [X for (X,Y) in brk[u,t]] bbrk = [Y for (X,Y) in brk[u,t]] # convex combination part: K = len(brk[u,t])-1 z = {} for k in range(K+1): z[k] = model.addVar(ub=1) # do not name variables for avoiding clash model.update() model.addConstr(p[u,t] == quicksum(abrk[k]*z[k] for k in range(K+1))) model.addConstr(fuel[u,t] == quicksum(bbrk[k]*z[k] for k in range(K+1))) model.addConstr(quicksum(z[k] for k in range(K+1)) == y[u,t]) model.addSOS(GRB.SOS_TYPE2, [z[k] for k in range(K+1)]) # objective model.addConstr(startup[u,t] == a_hot[u]*hot[u,t] + a_cold[u]*cold[u,t]) obj += fuel[u,t] + startup[u,t] model.setObjective(obj, GRB.MINIMIZE) def printsol(y, p, hot, cold, on, off, fuel, startup, objval): """printsol: display solution (for human consumption)""" print "Obj:", objval print "fuel:", sum(fuel[u,t].x for u in U for t in range(1,T+1)) print "startup:", sum(startup[u,t].x for u in U for t in range(1,T+1)) print "\ty, on, off:" for t in range(T+1): print t, "\t", for u in U: print int(y[u,t].x+.5), print "\t", if t == 0: print continue for u in U: print int(on[u,t].x+.5), print "\t", for u in U: print int(off[u,t].x+.5), print print "startup:" for t in range(1,T+1): print t,"\t", for u in U: print round(a_hot[u]*hot[u,t].x + a_cold[u]*cold[u,t].x,3), print "\t*", print round(sum(startup[u,t].x for u in U),3), print print "fuel:" for t in range(1,T+1): print t,"\t", for u in U: print round(fuel[u,t].x,3), print "\t*", print round(sum(fuel[u,t].x for u in U),3), print print "prod levels:" for t in range(1,T+1): print t, "\t", for u in U: print round(p[u,t].x,3), print "\t*", print round(sum(p[u,t].x for u in U),3), print if __name__ == "__main__": import sys, time try: cpu, inst, nperiods = float(sys.argv[1]), sys.argv[2], int(sys.argv[3]) except: print "usage: %s CPU inst nperiods" % sys.argv[0] print " where:" print " - CPU is the limit CPU allowed" print " - inst is one of: toy eld13 eld40 ucp10 ucp13 ucp40" print " - nperiods is the number of periods to consider, from 1 to 24" exit(0) # inst = "ucp10" # nperiods = 5 exec("U,T,a,b,c,e,f,p_min,p_max,y_prev,t_cold,n_init,min_on,min_off,a_hot,a_cold,dem,res = %s(%d)"%(inst,nperiods)) try: exec("U,T,a,b,c,e,f,p_min,p_max,y_prev,t_cold,n_init,min_on,min_off,a_hot,a_cold,dem,res = %s(%d)"%(inst,nperiods)) except: print "could not load instance %s(%d)"%(inst,nperiods) exit(0) brk = {} prev = {} # for u in U: for t in range(1,T+1): brk[u,t] = lower_brkpts(a[u],b[u],c[u],0,e[u],f[u],p_min[u],p_max[u],1) prev[u,t] = [] assert (p_max[u],cost(a[u],b[u],c[u],0,e[u],f[u],p_min[u],p_max[u])) in brk[u,t] EPS = 1.e-7 n = 100 has_solution = False strict = False while True: print "\nCURRENT n:", n lower = ucp_base(U, T, p_min, p_max, y_prev, t_cold, n_init, min_on, min_off, a_hot, a_cold, dem, res) y, p, hot, cold, on, off, fuel, startup = lower.__data brk = {} for u in U: print u,'len brk:', for t in range(1,T+1): brk[u,t] = expand_brkpts(a[u],b[u],c[u],0,e[u],f[u],p_min[u],p_max[u],prev[u,t],n) print len(brk[u,t]), print ucp_lin_obj(lower, U, T, brk) if has_solution: for u in U: for t in range(1,T+1): y[u,t].start = yp[u,t] hot[u,t].start = hotp[u,t] cold[u,t].start = coldp[u,t] on[u,t].start = onp[u,t] off[u,t].start = offp[u,t] lower.Params.Threads = 1 cpulim = cpu - time.clock() if cpulim < 0: print "time limit reached" break lower.Params.TimeLimit = cpulim if strict: lower.Params.MIPGap = 1.e-12 lower.Params.MIPGapAbs = 1.e-12 lower.Params.IntFeasTol = 1.e-9 lower.Params.FeasibilityTol = 1.e-9 # else: # does not work: default parameters are usually less strict than error/10 # try: # lower.Params.MIPGap = error/10 # except: # lower.Params.MIPGap = .1 lower.optimize() has_solution = True # printsol(y, p, hot, cold, on, off, fuel, startup, lower.ObjVal) fixed = sum(startup[u,t].x for u in U for t in range(1,T+1)) variable = sum(cost(a[u],b[u],c[u],0,e[u],f[u],p_min[u],p[u,t].x) for u in U for t in range(1,T+1) if y[u,t].x > .5) error = (fixed+variable - lower.ObjBound)/lower.ObjBound print "CURRENT ERROR: LB UP rel abs:", lower.ObjBound, "&", fixed+variable, print "&", error, "&", (fixed+variable - lower.ObjBound), "&", time.clock() if strict and error < EPS: break if error < EPS or error < lower.Params.MIPGap: strict = True if lower.status != GRB.Status.OPTIMAL: print "time limit reached" break yp, hotp, coldp, onp, offp = {}, {}, {}, {}, {} if has_solution: for u in U: for t in range(1,T+1): yp[u,t] = y[u,t].start hotp[u,t] = hot[u,t].start coldp[u,t] = cold[u,t].start onp[u,t] = on[u,t].start offp[u,t] = off[u,t].start # print "prev:", prev new = False for u in U: for t in range(1,T+1): if y[u,t].x > .5: # print u, t, p[u,t].x, "error:", cost(a[u],b[u],c[u],0,e[u],f[u],p_min[u],p[u,t].x), "-", fuel[u,t].x, \ # "=", cost(a[u],b[u],c[u],0,e[u],f[u],p_min[u],p[u,t].x) - fuel[u,t].x if has_added_p(p_min[u],p_max[u],f[u],prev[u,t],p[u,t].x): new = True # print "adding value for", u, t, p[u,t].x for pp in prev[u,t]: # print u, t, pp, p[u,t].x assert abs(pp-p[u,t].x) > 1.e-6 prev[u,t].append(p[u,t].x) prev[u,t].sort() else: # print "new breakpoint for", u, t, "not required" pass if not new: n *= 2 strict = True y, p, hot, cold, on, off, fuel, startup = lower.__data printsol(y, p, hot, cold, on, off, fuel, startup, lower.ObjVal)