Merge pull request #25 from jylambert/obj_vars

added variables to obj function, improved tests

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}')

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.")

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)

tests/mts_easy.py → 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

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

tests/test_setup.py

@@ -2,6 +2,7 @@
 
 from pytest import approx
 
+
 def test_import():
     import topotherm
 

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

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(