diff --git a/examples/one-producer/run_mts.py b/examples/one-producer/run_mts.py index a1762eb..d224d90 100644 --- a/examples/one-producer/run_mts.py +++ b/examples/one-producer/run_mts.py @@ -17,7 +17,7 @@ DATAPATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data_mts') OUTPUTPATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), - 'results/mts_eco_easy/') + 'results/mts_eco/') REGRESSION = 'regression.csv' # regression coefficients for thermal capacity and heat losses TIMESERIES = 'timeseries.csv' # timeseries for heat scaling PLOTS = True # save plots of the district @@ -110,6 +110,9 @@ def main(filepath, outputpath, plots=True, solver='gurobi', mode='forced'): if __name__ == '__main__': - main(filepath=os.path.join(DATAPATH), outputpath=os.path.join(OUTPUTPATH), - plots=PLOTS, solver=SOLVER, mode='economic') + main(filepath=os.path.join(DATAPATH), + outputpath=os.path.join(OUTPUTPATH), + plots=PLOTS, + solver=SOLVER, + mode='forced') print(f'Finished {OUTPUTPATH}') diff --git a/tests/conftest.py b/tests/conftest.py index ab4b8c1..62fcd5d 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,3 +1,9 @@ +"""Pytest configuration file.""" + + def pytest_addoption(parser): - parser.addoption("--solver", action="store", default="scip", - help="Define the solver to use for the optimization. Default is cbc.") + parser.addoption( + "--solver", + action="store", + default="scip", + help="Define the solver to use for the optimization. Default is cbc.") diff --git a/tests/hydraulics.py b/tests/hydraulics.py index 5a24976..53a01b0 100644 --- a/tests/hydraulics.py +++ b/tests/hydraulics.py @@ -34,6 +34,7 @@ def test_regression(): # heat loss regression r_heat_loss = precalc.regression_heat_losses( settings, thermal_capacity=r_thermal_cap) + # assert that each params item is within 1% of the expected value assert r_heat_loss['a'] == approx(4.348000e-07, rel=0.01) assert r_heat_loss['b'] == approx(0.02189, rel=0.01) diff --git a/tests/mts_easy.py b/tests/mts.py similarity index 90% rename from tests/mts_easy.py rename to tests/mts.py index b523626..942a6b1 100644 --- a/tests/mts_easy.py +++ b/tests/mts.py @@ -1,10 +1,14 @@ +"""Test the multiple timestep optimization model.""" + import os import pyomo.environ as pyo import pandas as pd +from pytest import approx import topotherm as tt + def read_regression(path, i): """Read the regression coefficients for the thermal capacity and heat losses from csv file. @@ -24,7 +28,7 @@ def read_regression(path, i): return r_thermal_cap, r_heat_loss -def test_mtseasy_forced(request): +def test_mts_forced(request): """Main function to run the optimization""" solver_name = request.config.getoption("--solver") assert solver_name in ["scip", "gurobi", "cbc"], f"Unsupported solver: {solver_name}" @@ -62,10 +66,10 @@ def test_mtseasy_forced(request): assert result.solver.status == pyo.SolverStatus.ok # assert that the objective value is less than 2% away from the expected # value - assert (abs(pyo.value(model.obj)) - 4.6259e+06) < 0.02 * 4.6259e+06 + assert abs(pyo.value(model.obj)) == approx(519676.4358105995, rel=0.02) -def test_mtseasy_eco(request): +def test_mts_eco(request): """Main function to run the optimization""" solver_name = request.config.getoption("--solver") assert solver_name in ["scip", "gurobi", "cbc"], f"Unsupported solver: {solver_name}" @@ -100,4 +104,5 @@ def test_mtseasy_eco(request): result = opt.solve(model, tee=True) assert result.solver.status == pyo.SolverStatus.ok - assert (abs(pyo.value(model.obj)) - 4.01854e+04) < 0.02 * 4.01854e+04 + assert pyo.value(model.obj) == approx(-39212.8646266408, rel=0.02) + # assert (abs(pyo.value(model.obj)) - 4.01854e+04) < 0.02 * 4.01854e+04 diff --git a/tests/sts.py b/tests/sts.py index df410b8..19fe5ea 100644 --- a/tests/sts.py +++ b/tests/sts.py @@ -62,7 +62,7 @@ def test_sts_forced(request): assert result.solver.status == pyo.SolverStatus.ok # assert that the objective value is less than 2% away from the expected # value - assert abs(pyo.value(model.obj)) == approx(519676.4358105995, rel=0.02) + assert pyo.value(model.obj) == approx(519676.4358105995, rel=0.02) def test_sts_eco(request): @@ -101,5 +101,5 @@ def test_sts_eco(request): result = opt.solve(model, tee=True) assert result.solver.status == pyo.SolverStatus.ok - assert abs(pyo.value(model.obj)) == approx(4.01854e+04, rel=0.02) + assert pyo.value(model.obj) == approx(-4.01854e+04, rel=0.02) # assert (abs(pyo.value(model.obj)) - 4.01854e+04) < 0.02 * 4.01854e+04 diff --git a/tests/test_setup.py b/tests/test_setup.py index fbe98ec..72d6c12 100644 --- a/tests/test_setup.py +++ b/tests/test_setup.py @@ -2,6 +2,7 @@ from pytest import approx + def test_import(): import topotherm diff --git a/topotherm/multiple_timestep.py b/topotherm/multiple_timestep.py index c207229..4f4b4c2 100644 --- a/topotherm/multiple_timestep.py +++ b/topotherm/multiple_timestep.py @@ -3,7 +3,8 @@ The module contains the following functions: * annuity: Calculate the annuity factor - * model: Create the optimization model for the multiple time steps operation + * model: Create the optimization model for the multiple time steps + operation """ import numpy as np import pyomo.environ as pyo @@ -260,56 +261,63 @@ def connection_to_consumer_fcd(m, d, j, t): rule=connection_to_consumer_fcd, doc=msg_) - if optimization_mode == "forced": def connection_to_consumer_built_fcd(m, j): return m.lambda_b[j] >= 1 mdl.cons_connection_forced = pyo.Constraint( - mdl.forced_edges, rule=connection_to_consumer_built_fcd, + mdl.forced_edges, + rule=connection_to_consumer_built_fcd, doc='The house connection has to be built to satisfy the demand') - if optimization_mode == "forced": def total_energy_conservation(m, t): return sum(m.P_source[k, t] for k in m.set_n_p) \ - sum(m.P['ij', 'in', k, t] - m.P['ij', 'out', k, t] for k in m.set_n_i) \ - sum(m.P['ji', 'in', k, t] - m.P['ji', 'out', k, t] for k in m.set_n_i) \ - sum(matrices['q_c'][k, t] for k in m.set_n_c) == 0 - mdl.cons_total_energy_cons = pyo.Constraint(mdl.set_t, rule=total_energy_conservation, - doc='Total energy conservation') - - def objective_function(m): - fuel = sum( - sum(m.P_source[k, t] + mdl.cons_total_energy_cons = pyo.Constraint( + mdl.set_t, + rule=total_energy_conservation, + doc='Total energy conservation') + + mdl.revenue = pyo.Var(doc='Revenue', domain=pyo.NegativeReals) + mdl.revenue_constr = pyo.Constraint( + expr=mdl.revenue == sum( + sum(mdl.lambda_b[j] * matrices['flh_consumer'][k, t] + * matrices['q_c'][k, t] for k, j in mdl.consumer_edges) + for t in mdl.set_t) * economics.heat_price * (-1), + doc='Revenue constraint') + + mdl.opex_source = pyo.Var(doc='OPEX Source', domain=pyo.PositiveReals) + mdl.opex_source_constr = pyo.Constraint( + expr=mdl.opex_source == sum( + sum(mdl.P_source[k, t] * economics.source_price[k] * matrices['flh_source'][k, t] - for k in m.set_n_p) - for t in mdl.set_t) - - pipes = sum( - ( - (m.P_cap[k] * regression_inst['a'] - + regression_inst['b'] * m.lambda_b[k] - ) * annuity(economics.pipes_c_irr, economics.pipes_lifetime) * matrices['l_i'][k] - ) for k in m.set_n_i - ) - - source = sum(m.P_source_inst[k] + for k in mdl.set_n_p) + for t in mdl.set_t), + doc='OPEX Source constraint') + + mdl.capex_pipes = pyo.Var(doc='CAPEX Pipe', domain=pyo.PositiveReals) + pipes = sum( + ((mdl.P_cap[k] * regression_inst['a'] + + regression_inst['b'] * mdl.lambda_b[k]) + * annuity(economics.pipes_c_irr, economics.pipes_lifetime) + * matrices['l_i'][k]) for k in mdl.set_n_i) + + mdl.capex_pipe_constr = pyo.Constraint( + expr=mdl.capex_pipes == pipes, + doc='CAPEX Pipe constraint') + + mdl.capex_source = pyo.Var(doc='CAPEX Source', domain=pyo.PositiveReals) + mdl.capex_source_constr = pyo.Constraint( + expr=mdl.capex_source == sum(mdl.P_source_inst[k] * economics.source_c_inv[k] * annuity(economics.source_c_irr[k], economics.source_lifetime[k]) - for k in m.set_n_p) - - if optimization_mode == "economic": - revenue = sum( - sum(m.lambda_b[j] - * matrices['flh_consumer'][k, t] - * matrices['q_c'][k, t] - for k, j in mdl.consumer_edges) - for t in mdl.set_t) * economics.heat_price * (-1) - else: - revenue = 0 - - return fuel + pipes + source + revenue + for k in mdl.set_n_p), + doc='CAPEX Source constraint') - mdl.obj = pyo.Objective(rule=objective_function, doc='Objective function') + mdl.obj = pyo.Objective( + expr=mdl.capex_source + mdl.capex_pipes + mdl.opex_source + mdl.revenue, + doc='Objective function') return mdl diff --git a/topotherm/single_timestep.py b/topotherm/single_timestep.py index 2f42100..8a60589 100644 --- a/topotherm/single_timestep.py +++ b/topotherm/single_timestep.py @@ -247,7 +247,7 @@ def total_energy_cons(m, t): mdl.set_t, rule=total_energy_cons, doc='Total energy conservation') - mdl.revenue = pyo.Var(doc='Revenue') + mdl.revenue = pyo.Var(doc='Revenue', domain=pyo.NegativeReals) mdl.revenue_constr = pyo.Constraint( expr=mdl.revenue == ( sum( @@ -267,7 +267,7 @@ def total_energy_cons(m, t): for t in mdl.set_t) * economics.heat_price * (-1)), doc='Revenue constraint') - mdl.opex_source = pyo.Var(doc='OPEX Source') + mdl.opex_source = pyo.Var(doc='OPEX Source', domain=pyo.PositiveReals) mdl.opex_source_constr = pyo.Constraint( expr=mdl.opex_source == sum( sum(mdl.P_source[k, t] @@ -277,7 +277,7 @@ def total_energy_cons(m, t): for t in mdl.set_t), doc='OPEX Source constraint') - mdl.capex_pipes = pyo.Var(doc='CAPEX Pipe') + mdl.capex_pipes = pyo.Var(doc='CAPEX Pipe', domain=pyo.PositiveReals) # CAREFUL HARDCODED FOR 0 TIME STEPS def pipes_fix(k): @@ -289,19 +289,18 @@ def pipes_var(k): * (mdl.lambda_['ij', k] + mdl.lambda_['ji', k])) mdl.capex_pipe_constr = pyo.Constraint( - expr=mdl.capex_pipes == (sum( + expr=mdl.capex_pipes == sum( (pipes_fix(k) + pipes_var(k)) * matrices['l_i'][k] for k in mdl.set_n_i) - * annuity(economics.pipes_c_irr, economics.pipes_lifetime)), + * annuity(economics.pipes_c_irr, economics.pipes_lifetime), doc='CAPEX Pipe constraint') - mdl.capex_source = pyo.Var(doc='CAPEX Source') + mdl.capex_source = pyo.Var(doc='CAPEX Source', domain=pyo.PositiveReals) mdl.capex_source_constr = pyo.Constraint( - expr=mdl.capex_source == sum(mdl.P_source_inst[k] - * economics.source_c_inv[k] - * annuity(economics.source_c_irr[k], - economics.source_lifetime[k]) - for k in mdl.set_n_p), + expr=mdl.capex_source == sum( + mdl.P_source_inst[k] * economics.source_c_inv[k] + * annuity(economics.source_c_irr[k], economics.source_lifetime[k]) + for k in mdl.set_n_p), doc='CAPEX Source constraint') mdl.obj = pyo.Objective(