|
| 1 | +from pyomo.environ import * |
| 2 | +import json |
| 3 | +import sys |
| 4 | + |
| 5 | +## Grab instance file from first command line argument |
| 6 | +data_file = sys.argv[1] |
| 7 | + |
| 8 | +print('loading data') |
| 9 | +data = json.load(open(data_file, 'r')) |
| 10 | + |
| 11 | +thermal_gens = data['thermal_generators'] |
| 12 | +renewable_gens = data['renewable_generators'] |
| 13 | + |
| 14 | +time_periods = {t+1 : t for t in range(data['time_periods'])} |
| 15 | + |
| 16 | +gen_startup_categories = {g : list(range(0, len(gen['startup']))) for (g, gen) in thermal_gens.items()} |
| 17 | +gen_pwl_points = {g : list(range(0, len(gen['piecewise_production']))) for (g, gen) in thermal_gens.items()} |
| 18 | + |
| 19 | +print('building model') |
| 20 | +m = ConcreteModel() |
| 21 | + |
| 22 | +m.cg = Var(thermal_gens.keys(), time_periods.keys()) |
| 23 | +m.pg = Var(thermal_gens.keys(), time_periods.keys(), within=NonNegativeReals) |
| 24 | +m.rg = Var(thermal_gens.keys(), time_periods.keys(), within=NonNegativeReals) |
| 25 | +m.pw = Var(renewable_gens.keys(), time_periods.keys(), within=NonNegativeReals) |
| 26 | +m.ug = Var(thermal_gens.keys(), time_periods.keys(), within=Binary) |
| 27 | +m.vg = Var(thermal_gens.keys(), time_periods.keys(), within=Binary) |
| 28 | +m.wg = Var(thermal_gens.keys(), time_periods.keys(), within=Binary) |
| 29 | + |
| 30 | +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) ## |
| 31 | +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) ## |
| 32 | + |
| 33 | +m.obj = Objective(expr=sum( |
| 34 | + sum( |
| 35 | + m.cg[g,t] + gen['piecewise_production'][0]['cost']*m.ug[g,t] |
| 36 | + + sum( gen_startup['cost']*m.dg[g,s,t] for (s, gen_startup) in enumerate(gen['startup'])) |
| 37 | + for t in time_periods) |
| 38 | + for g, gen in thermal_gens.items() ) |
| 39 | + ) #(1) |
| 40 | + |
| 41 | +m.demand = Constraint(time_periods.keys()) |
| 42 | +m.reserves = Constraint(time_periods.keys()) |
| 43 | +for t,t_idx in time_periods.items(): |
| 44 | + 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) |
| 45 | + m.reserves[t] = sum( m.rg[g,t] for g in thermal_gens ) >= data['reserves'][t_idx] #(3) |
| 46 | + |
| 47 | +m.uptimet0 = Constraint(thermal_gens.keys()) |
| 48 | +m.downtimet0 = Constraint(thermal_gens.keys()) |
| 49 | +m.logicalt0 = Constraint(thermal_gens.keys()) |
| 50 | +m.startupt0 = Constraint(thermal_gens.keys()) |
| 51 | + |
| 52 | +m.rampupt0 = Constraint(thermal_gens.keys()) |
| 53 | +m.rampdownt0 = Constraint(thermal_gens.keys()) |
| 54 | +m.shutdownt0 = Constraint(thermal_gens.keys()) |
| 55 | + |
| 56 | +for g, gen in thermal_gens.items(): |
| 57 | + if gen['unit_on_t0'] == 1: |
| 58 | + if gen['time_up_minimum'] - gen['time_up_t0'] >= 1: |
| 59 | + 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) |
| 60 | + elif gen['unit_on_t0'] == 0: |
| 61 | + if gen['time_down_minimum'] - gen['time_down_t0'] >= 1: |
| 62 | + 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) |
| 63 | + else: |
| 64 | + raise Exception('Invalid unit_on_t0 for generator {}, unit_on_t0={}'.format(g, gen['unit_on_t0'])) |
| 65 | + |
| 66 | + m.logicalt0[g] = m.ug[g,1] - gen['unit_on_t0'] == m.vg[g,1] - m.wg[g,1] #(6) |
| 67 | + |
| 68 | + startup_expr = sum( |
| 69 | + sum( m.dg[g,s,t] |
| 70 | + for t in range( |
| 71 | + max(1,gen['startup'][s+1]['lag']-gen['time_down_t0']+1), |
| 72 | + min(gen['startup'][s+1]['lag']-1,data['time_periods'])+1 |
| 73 | + ) |
| 74 | + ) |
| 75 | + for s,_ in enumerate(gen['startup'][:-1])) ## all but last |
| 76 | + if isinstance(startup_expr, int): |
| 77 | + pass |
| 78 | + else: |
| 79 | + m.startupt0[g] = startup_expr == 0 #(7) |
| 80 | + |
| 81 | + 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) |
| 82 | + |
| 83 | + m.rampdownt0[g] = gen['unit_on_t0']*(gen['power_output_t0']-gen['power_output_minimum']) - m.pg[g,1] <= gen['ramp_down_limit'] #(9) |
| 84 | + |
| 85 | + |
| 86 | + 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) |
| 87 | + |
| 88 | + if isinstance(shutdown_constr, bool): |
| 89 | + pass |
| 90 | + else: |
| 91 | + m.shutdownt0[g] = shutdown_constr |
| 92 | + |
| 93 | +m.mustrun = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 94 | +m.logical = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 95 | +m.uptime = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 96 | +m.downtime = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 97 | +m.startup_select = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 98 | +m.gen_limit1 = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 99 | +m.gen_limit2 = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 100 | +m.ramp_up = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 101 | +m.ramp_down = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 102 | +m.power_select = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 103 | +m.cost_select = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 104 | +m.on_select = Constraint(thermal_gens.keys(), time_periods.keys()) |
| 105 | + |
| 106 | +for g, gen in thermal_gens.items(): |
| 107 | + for t in time_periods: |
| 108 | + m.mustrun[g,t] = m.ug[g,t] >= gen['must_run'] #(11) |
| 109 | + |
| 110 | + if t > 1: |
| 111 | + m.logical[g,t] = m.ug[g,t] - m.ug[g,t-1] == m.vg[g,t] - m.wg[g,t] #(12) |
| 112 | + |
| 113 | + UT = min(gen['time_up_minimum'],data['time_periods']) |
| 114 | + if t >= UT: |
| 115 | + m.uptime[g,t] = sum(m.vg[g,t] for t in range(t-UT+1, t+1)) <= m.ug[g,t] #(13) |
| 116 | + DT = min(gen['time_down_minimum'],data['time_periods']) |
| 117 | + if t >= DT: |
| 118 | + m.downtime[g,t] = sum(m.wg[g,t] for t in range(t-DT+1, t+1)) <= 1-m.ug[g,t] #(14) |
| 119 | + m.startup_select[g,t] = m.vg[g,t] == sum(m.dg[g,s,t] for s,_ in enumerate(gen['startup'])) #(16) |
| 120 | + |
| 121 | + 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) |
| 122 | + |
| 123 | + if t < len(time_periods): |
| 124 | + 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) |
| 125 | + |
| 126 | + if t > 1: |
| 127 | + m.ramp_up[g,t] = m.pg[g,t]+m.rg[g,t] - m.pg[g,t-1] <= gen['ramp_up_limit'] #(19) |
| 128 | + m.ramp_down[g,t] = m.pg[g,t-1] - m.pg[g,t] <= gen['ramp_down_limit'] #(20 |
| 129 | + |
| 130 | + piece_mw1 = gen['piecewise_production'][0]['mw'] |
| 131 | + piece_cost1 = gen['piecewise_production'][0]['cost'] |
| 132 | + 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) |
| 133 | + 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) |
| 134 | + m.on_select[g,t] = m.ug[g,t] == sum(m.lg[g,l,t] for l,_ in enumerate(gen['piecewise_production'])) #(23) |
| 135 | + |
| 136 | +m.startup_allowed = Constraint(m.dg_index) |
| 137 | +for g, gen in thermal_gens.items(): |
| 138 | + for s,_ in enumerate(gen['startup'][:-1]): ## all but last |
| 139 | + for t in time_periods: |
| 140 | + if t >= gen['startup'][s+1]['lag']: |
| 141 | + 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) |
| 142 | + |
| 143 | +for w, gen in renewable_gens.items(): |
| 144 | + for t, t_idx in time_periods.items(): |
| 145 | + m.pw[w,t].setlb(gen['power_output_minimum'][t_idx]) #(24) |
| 146 | + m.pw[w,t].setub(gen['power_output_maximum'][t_idx]) #(24) |
| 147 | + |
| 148 | +print("model setup complete") |
| 149 | + |
| 150 | +from pyomo.opt import SolverFactory |
| 151 | +cbc = SolverFactory('cbc') |
| 152 | + |
| 153 | +print("solving") |
| 154 | +cbc.solve(m, options={'ratioGap':0.01}, tee=True) |
0 commit comments