From 2cb5e08c8f0c50fc0bcadade05affef8de517d58 Mon Sep 17 00:00:00 2001 From: whart222 Date: Mon, 10 Jun 2024 09:12:05 -0600 Subject: [PATCH 1/3] Several changes 1. Fixed error when referencing an index set 2. Enabled command-line specification of solver --- uc_model.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/uc_model.py b/uc_model.py index 028e640..7030386 100644 --- a/uc_model.py +++ b/uc_model.py @@ -3,7 +3,14 @@ import sys ## Grab instance file from first command line argument +if len(sys.argv) == 1: + print("uc_model.py []") + sys.exit(0) data_file = sys.argv[1] +if len(sys.argv) == 3: + solver_name = sys.argv[2] +else: + solver_name = 'cbc' print('loading data') data = json.load(open(data_file, 'r')) @@ -133,7 +140,7 @@ m.cost_select[g,t] = m.cg[g,t] == sum( (piece['cost'] - piece_cost1)*m.lg[g,l,t] for l,piece in enumerate(gen['piecewise_production'])) #(22) m.on_select[g,t] = m.ug[g,t] == sum(m.lg[g,l,t] for l,_ in enumerate(gen['piecewise_production'])) #(23) -m.startup_allowed = Constraint(m.dg_index) +m.startup_allowed = Constraint(m.dg.index_set()) for g, gen in thermal_gens.items(): for s,_ in enumerate(gen['startup'][:-1]): ## all but last for t in time_periods: @@ -148,7 +155,9 @@ print("model setup complete") from pyomo.opt import SolverFactory -cbc = SolverFactory('cbc') +solver = SolverFactory(solver_name) print("solving") -cbc.solve(m, options={'ratioGap':0.01}, tee=True) +options = {'gurobi':{'ratioGap':0.01}, 'glpk':{'mipgap':0.01}, 'cbc':{'ratioGap':0.01}, 'scip':{'mipgap':0.01}} +solver.solve(m, options=options[solver_name], tee=True) + From 5f829ac552163a4fb3d56afa72cf68b5b31de4f4 Mon Sep 17 00:00:00 2001 From: whart222 Date: Mon, 10 Jun 2024 09:13:34 -0600 Subject: [PATCH 2/3] Reformatting with black --- uc_model.py | 281 +++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 200 insertions(+), 81 deletions(-) diff --git a/uc_model.py b/uc_model.py index 7030386..d0c7be8 100644 --- a/uc_model.py +++ b/uc_model.py @@ -10,46 +10,78 @@ if len(sys.argv) == 3: solver_name = sys.argv[2] else: - solver_name = 'cbc' + solver_name = "cbc" -print('loading data') -data = json.load(open(data_file, 'r')) +print("loading data") +data = json.load(open(data_file, "r")) -thermal_gens = data['thermal_generators'] -renewable_gens = data['renewable_generators'] +thermal_gens = data["thermal_generators"] +renewable_gens = data["renewable_generators"] -time_periods = {t+1 : t for t in range(data['time_periods'])} +time_periods = {t + 1: t for t in range(data["time_periods"])} -gen_startup_categories = {g : list(range(0, len(gen['startup']))) for (g, gen) in thermal_gens.items()} -gen_pwl_points = {g : list(range(0, len(gen['piecewise_production']))) for (g, gen) in thermal_gens.items()} +gen_startup_categories = { + g: list(range(0, len(gen["startup"]))) for (g, gen) in thermal_gens.items() +} +gen_pwl_points = { + g: list(range(0, len(gen["piecewise_production"]))) + for (g, gen) in thermal_gens.items() +} -print('building model') +print("building model") m = ConcreteModel() m.cg = Var(thermal_gens.keys(), time_periods.keys()) -m.pg = Var(thermal_gens.keys(), time_periods.keys(), within=NonNegativeReals) -m.rg = Var(thermal_gens.keys(), time_periods.keys(), within=NonNegativeReals) +m.pg = Var(thermal_gens.keys(), time_periods.keys(), within=NonNegativeReals) +m.rg = Var(thermal_gens.keys(), time_periods.keys(), within=NonNegativeReals) m.pw = Var(renewable_gens.keys(), time_periods.keys(), within=NonNegativeReals) -m.ug = Var(thermal_gens.keys(), time_periods.keys(), within=Binary) -m.vg = Var(thermal_gens.keys(), time_periods.keys(), within=Binary) -m.wg = Var(thermal_gens.keys(), time_periods.keys(), within=Binary) - -m.dg = Var(((g,s,t) for g in thermal_gens for s in gen_startup_categories[g] for t in time_periods), within=Binary) ## -m.lg = Var(((g,l,t) for g in thermal_gens for l in gen_pwl_points[g] for t in time_periods), within=UnitInterval) ## - -m.obj = Objective(expr=sum( - sum( - m.cg[g,t] + gen['piecewise_production'][0]['cost']*m.ug[g,t] - + sum( gen_startup['cost']*m.dg[g,s,t] for (s, gen_startup) in enumerate(gen['startup'])) - for t in time_periods) - for g, gen in thermal_gens.items() ) - ) #(1) +m.ug = Var(thermal_gens.keys(), time_periods.keys(), within=Binary) +m.vg = Var(thermal_gens.keys(), time_periods.keys(), within=Binary) +m.wg = Var(thermal_gens.keys(), time_periods.keys(), within=Binary) + +m.dg = Var( + ( + (g, s, t) + for g in thermal_gens + for s in gen_startup_categories[g] + for t in time_periods + ), + within=Binary, +) ## +m.lg = Var( + ((g, l, t) for g in thermal_gens for l in gen_pwl_points[g] for t in time_periods), + within=UnitInterval, +) ## + +m.obj = Objective( + expr=sum( + sum( + m.cg[g, t] + + gen["piecewise_production"][0]["cost"] * m.ug[g, t] + + sum( + gen_startup["cost"] * m.dg[g, s, t] + for (s, gen_startup) in enumerate(gen["startup"]) + ) + for t in time_periods + ) + for g, gen in thermal_gens.items() + ) +) # (1) m.demand = Constraint(time_periods.keys()) m.reserves = Constraint(time_periods.keys()) -for t,t_idx in time_periods.items(): - m.demand[t] = sum( m.pg[g,t]+gen['power_output_minimum']*m.ug[g,t] for (g, gen) in thermal_gens.items() ) + sum( m.pw[w,t] for w in renewable_gens ) == data['demand'][t_idx] #(2) - m.reserves[t] = sum( m.rg[g,t] for g in thermal_gens ) >= data['reserves'][t_idx] #(3) +for t, t_idx in time_periods.items(): + m.demand[t] = ( + sum( + m.pg[g, t] + gen["power_output_minimum"] * m.ug[g, t] + for (g, gen) in thermal_gens.items() + ) + + sum(m.pw[w, t] for w in renewable_gens) + == data["demand"][t_idx] + ) # (2) + m.reserves[t] = ( + sum(m.rg[g, t] for g in thermal_gens) >= data["reserves"][t_idx] + ) # (3) m.uptimet0 = Constraint(thermal_gens.keys()) m.downtimet0 = Constraint(thermal_gens.keys()) @@ -61,36 +93,82 @@ m.shutdownt0 = Constraint(thermal_gens.keys()) for g, gen in thermal_gens.items(): - if gen['unit_on_t0'] == 1: - if gen['time_up_minimum'] - gen['time_up_t0'] >= 1: - m.uptimet0[g] = sum( (m.ug[g,t] - 1) for t in range(1, min(gen['time_up_minimum'] - gen['time_up_t0'], data['time_periods'])+1)) == 0 #(4) - elif gen['unit_on_t0'] == 0: - if gen['time_down_minimum'] - gen['time_down_t0'] >= 1: - m.downtimet0[g] = sum( m.ug[g,t] for t in range(1, min(gen['time_down_minimum'] - gen['time_down_t0'], data['time_periods'])+1)) == 0 #(5) + if gen["unit_on_t0"] == 1: + if gen["time_up_minimum"] - gen["time_up_t0"] >= 1: + m.uptimet0[g] = ( + sum( + (m.ug[g, t] - 1) + for t in range( + 1, + min( + gen["time_up_minimum"] - gen["time_up_t0"], + data["time_periods"], + ) + + 1, + ) + ) + == 0 + ) # (4) + elif gen["unit_on_t0"] == 0: + if gen["time_down_minimum"] - gen["time_down_t0"] >= 1: + m.downtimet0[g] = ( + sum( + m.ug[g, t] + for t in range( + 1, + min( + gen["time_down_minimum"] - gen["time_down_t0"], + data["time_periods"], + ) + + 1, + ) + ) + == 0 + ) # (5) else: - raise Exception('Invalid unit_on_t0 for generator {}, unit_on_t0={}'.format(g, gen['unit_on_t0'])) - - m.logicalt0[g] = m.ug[g,1] - gen['unit_on_t0'] == m.vg[g,1] - m.wg[g,1] #(6) - - startup_expr = sum( - sum( m.dg[g,s,t] - for t in range( - max(1,gen['startup'][s+1]['lag']-gen['time_down_t0']+1), - min(gen['startup'][s+1]['lag']-1,data['time_periods'])+1 - ) - ) - for s,_ in enumerate(gen['startup'][:-1])) ## all but last + raise Exception( + "Invalid unit_on_t0 for generator {}, unit_on_t0={}".format( + g, gen["unit_on_t0"] + ) + ) + + m.logicalt0[g] = m.ug[g, 1] - gen["unit_on_t0"] == m.vg[g, 1] - m.wg[g, 1] # (6) + + startup_expr = sum( + sum( + m.dg[g, s, t] + for t in range( + max(1, gen["startup"][s + 1]["lag"] - gen["time_down_t0"] + 1), + min(gen["startup"][s + 1]["lag"] - 1, data["time_periods"]) + 1, + ) + ) + for s, _ in enumerate(gen["startup"][:-1]) + ) ## all but last if isinstance(startup_expr, int): pass else: - m.startupt0[g] = startup_expr == 0 #(7) - - m.rampupt0[g] = m.pg[g,1] + m.rg[g,1] - gen['unit_on_t0']*(gen['power_output_t0']-gen['power_output_minimum']) <= gen['ramp_up_limit'] #(8) - - m.rampdownt0[g] = gen['unit_on_t0']*(gen['power_output_t0']-gen['power_output_minimum']) - m.pg[g,1] <= gen['ramp_down_limit'] #(9) - - - shutdown_constr = gen['unit_on_t0']*(gen['power_output_t0']-gen['power_output_minimum']) <= gen['unit_on_t0']*(gen['power_output_maximum'] - gen['power_output_minimum']) - max((gen['power_output_maximum'] - gen['ramp_shutdown_limit']),0)*m.wg[g,1] #(10) + m.startupt0[g] = startup_expr == 0 # (7) + + m.rampupt0[g] = ( + m.pg[g, 1] + + m.rg[g, 1] + - gen["unit_on_t0"] * (gen["power_output_t0"] - gen["power_output_minimum"]) + <= gen["ramp_up_limit"] + ) # (8) + + m.rampdownt0[g] = ( + gen["unit_on_t0"] * (gen["power_output_t0"] - gen["power_output_minimum"]) + - m.pg[g, 1] + <= gen["ramp_down_limit"] + ) # (9) + + shutdown_constr = ( + gen["unit_on_t0"] * (gen["power_output_t0"] - gen["power_output_minimum"]) + <= gen["unit_on_t0"] + * (gen["power_output_maximum"] - gen["power_output_minimum"]) + - max((gen["power_output_maximum"] - gen["ramp_shutdown_limit"]), 0) + * m.wg[g, 1] + ) # (10) if isinstance(shutdown_constr, bool): pass @@ -112,52 +190,93 @@ for g, gen in thermal_gens.items(): for t in time_periods: - m.mustrun[g,t] = m.ug[g,t] >= gen['must_run'] #(11) + m.mustrun[g, t] = m.ug[g, t] >= gen["must_run"] # (11) if t > 1: - m.logical[g,t] = m.ug[g,t] - m.ug[g,t-1] == m.vg[g,t] - m.wg[g,t] #(12) + m.logical[g, t] = ( + m.ug[g, t] - m.ug[g, t - 1] == m.vg[g, t] - m.wg[g, t] + ) # (12) - UT = min(gen['time_up_minimum'],data['time_periods']) + UT = min(gen["time_up_minimum"], data["time_periods"]) if t >= UT: - m.uptime[g,t] = sum(m.vg[g,t] for t in range(t-UT+1, t+1)) <= m.ug[g,t] #(13) - DT = min(gen['time_down_minimum'],data['time_periods']) + m.uptime[g, t] = ( + sum(m.vg[g, t] for t in range(t - UT + 1, t + 1)) <= m.ug[g, t] + ) # (13) + DT = min(gen["time_down_minimum"], data["time_periods"]) if t >= DT: - m.downtime[g,t] = sum(m.wg[g,t] for t in range(t-DT+1, t+1)) <= 1-m.ug[g,t] #(14) - m.startup_select[g,t] = m.vg[g,t] == sum(m.dg[g,s,t] for s,_ in enumerate(gen['startup'])) #(16) - - m.gen_limit1[g,t] = m.pg[g,t]+m.rg[g,t] <= (gen['power_output_maximum'] - gen['power_output_minimum'])*m.ug[g,t] - max((gen['power_output_maximum'] - gen['ramp_startup_limit']),0)*m.vg[g,t] #(17) - - if t < len(time_periods): - m.gen_limit2[g,t] = m.pg[g,t]+m.rg[g,t] <= (gen['power_output_maximum'] - gen['power_output_minimum'])*m.ug[g,t] - max((gen['power_output_maximum'] - gen['ramp_shutdown_limit']),0)*m.wg[g,t+1] #(18) + m.downtime[g, t] = ( + sum(m.wg[g, t] for t in range(t - DT + 1, t + 1)) <= 1 - m.ug[g, t] + ) # (14) + m.startup_select[g, t] = m.vg[g, t] == sum( + m.dg[g, s, t] for s, _ in enumerate(gen["startup"]) + ) # (16) + + m.gen_limit1[g, t] = ( + m.pg[g, t] + m.rg[g, t] + <= (gen["power_output_maximum"] - gen["power_output_minimum"]) * m.ug[g, t] + - max((gen["power_output_maximum"] - gen["ramp_startup_limit"]), 0) + * m.vg[g, t] + ) # (17) + + if t < len(time_periods): + m.gen_limit2[g, t] = ( + m.pg[g, t] + m.rg[g, t] + <= (gen["power_output_maximum"] - gen["power_output_minimum"]) + * m.ug[g, t] + - max((gen["power_output_maximum"] - gen["ramp_shutdown_limit"]), 0) + * m.wg[g, t + 1] + ) # (18) if t > 1: - m.ramp_up[g,t] = m.pg[g,t]+m.rg[g,t] - m.pg[g,t-1] <= gen['ramp_up_limit'] #(19) - m.ramp_down[g,t] = m.pg[g,t-1] - m.pg[g,t] <= gen['ramp_down_limit'] #(20 - - piece_mw1 = gen['piecewise_production'][0]['mw'] - piece_cost1 = gen['piecewise_production'][0]['cost'] - m.power_select[g,t] = m.pg[g,t] == sum( (piece['mw'] - piece_mw1)*m.lg[g,l,t] for l,piece in enumerate(gen['piecewise_production'])) #(21) - m.cost_select[g,t] = m.cg[g,t] == sum( (piece['cost'] - piece_cost1)*m.lg[g,l,t] for l,piece in enumerate(gen['piecewise_production'])) #(22) - m.on_select[g,t] = m.ug[g,t] == sum(m.lg[g,l,t] for l,_ in enumerate(gen['piecewise_production'])) #(23) + m.ramp_up[g, t] = ( + m.pg[g, t] + m.rg[g, t] - m.pg[g, t - 1] <= gen["ramp_up_limit"] + ) # (19) + m.ramp_down[g, t] = ( + m.pg[g, t - 1] - m.pg[g, t] <= gen["ramp_down_limit"] + ) # (20 + + piece_mw1 = gen["piecewise_production"][0]["mw"] + piece_cost1 = gen["piecewise_production"][0]["cost"] + m.power_select[g, t] = m.pg[g, t] == sum( + (piece["mw"] - piece_mw1) * m.lg[g, l, t] + for l, piece in enumerate(gen["piecewise_production"]) + ) # (21) + m.cost_select[g, t] = m.cg[g, t] == sum( + (piece["cost"] - piece_cost1) * m.lg[g, l, t] + for l, piece in enumerate(gen["piecewise_production"]) + ) # (22) + m.on_select[g, t] = m.ug[g, t] == sum( + m.lg[g, l, t] for l, _ in enumerate(gen["piecewise_production"]) + ) # (23) m.startup_allowed = Constraint(m.dg.index_set()) for g, gen in thermal_gens.items(): - for s,_ in enumerate(gen['startup'][:-1]): ## all but last + for s, _ in enumerate(gen["startup"][:-1]): ## all but last for t in time_periods: - if t >= gen['startup'][s+1]['lag']: - m.startup_allowed[g,s,t] = m.dg[g,s,t] <= sum(m.wg[g,t-i] for i in range(gen['startup'][s]['lag'], gen['startup'][s+1]['lag'])) #(15) + if t >= gen["startup"][s + 1]["lag"]: + m.startup_allowed[g, s, t] = m.dg[g, s, t] <= sum( + m.wg[g, t - i] + for i in range( + gen["startup"][s]["lag"], gen["startup"][s + 1]["lag"] + ) + ) # (15) for w, gen in renewable_gens.items(): for t, t_idx in time_periods.items(): - m.pw[w,t].setlb(gen['power_output_minimum'][t_idx]) #(24) - m.pw[w,t].setub(gen['power_output_maximum'][t_idx]) #(24) + m.pw[w, t].setlb(gen["power_output_minimum"][t_idx]) # (24) + m.pw[w, t].setub(gen["power_output_maximum"][t_idx]) # (24) print("model setup complete") from pyomo.opt import SolverFactory + solver = SolverFactory(solver_name) print("solving") -options = {'gurobi':{'ratioGap':0.01}, 'glpk':{'mipgap':0.01}, 'cbc':{'ratioGap':0.01}, 'scip':{'mipgap':0.01}} +options = { + "gurobi": {"ratioGap": 0.01}, + "glpk": {"mipgap": 0.01}, + "cbc": {"ratioGap": 0.01}, + "scip": {"mipgap": 0.01}, +} solver.solve(m, options=options[solver_name], tee=True) - From a428814e3144f8980622a704c296dd9a7aa324d2 Mon Sep 17 00:00:00 2001 From: whart222 Date: Tue, 11 Jun 2024 05:36:27 -0600 Subject: [PATCH 3/3] Updating gurobi mipgap parameter --- uc_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uc_model.py b/uc_model.py index d0c7be8..f82fe55 100644 --- a/uc_model.py +++ b/uc_model.py @@ -274,7 +274,7 @@ print("solving") options = { - "gurobi": {"ratioGap": 0.01}, + "gurobi": {"mipgap": 0.01}, "glpk": {"mipgap": 0.01}, "cbc": {"ratioGap": 0.01}, "scip": {"mipgap": 0.01},