Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement alternative back-end to Pyomo (e.g. linopy) #1064

Open
fwitte opened this issue Mar 26, 2024 · 9 comments
Open

Implement alternative back-end to Pyomo (e.g. linopy) #1064

fwitte opened this issue Mar 26, 2024 · 9 comments

Comments

@fwitte
Copy link
Member

fwitte commented Mar 26, 2024

Or, potentially even replace pyomo (also see https://forum.openmod.org/t/linopy-for-oemof/4638).

I have recently used oemof-solph in a small dispatch optimization problem featuring two Converters, a Storage and a GenericCHP. I profiled the execution times and roughly found:

  • Model building: 5 s
  • Solving (with gurobi): 15 s
  • Results processing. 5 s
  • Some more pyomo and oemof-solph related overhead: 5 s

(not the actual times, the point is, that only half of the time is spent on solving)

There are a couple of issues with this for my use case:

  • About the solving, I obviously cannot do anything except maybe removing integer variables or clustering the time series (I do not necessarily want to do that at this moment).
  • I needed to modify the storage capacity and one of the converter nominal power values a couple of times, therefore needed to rebuild the model every time I wanted to solve.
  • I only wanted to have one time series from the results, it would be nice, if that could be extracted more time efficient.

Since I loose so much time in tasks which is not solving, I would like to bring this topic to the table.

I tested a very simple model (Source -> Bus -> Converter -> Bus -> Sink) in oemof-solph and linopy and got the following results (based on a single execution, I will retry a couple of times, so maybe do not look at the absolutes but more at the relative differences):

oemof-solph

from oemof import solph

es = solph.EnergySystem(timeindex=solph.create_time_index(2012, number=8761))

b1 = solph.Bus(label="b1")
b2 = solph.Bus(label="b2")

source = solph.components.Source(
    label="source", outputs={b1: solph.Flow()}
)
sink = solph.components.Sink(
    label="sink", inputs={b2: solph.Flow(nominal_value=1, fix=1)}
)

converter = solph.components.Converter(
    label="converter",
    inputs={b1: solph.Flow()},
    outputs={b2: solph.Flow(nominal_value=10)},
    conversion_factors={b2: 1}
)

es.add(b1, b2, source, sink, converter)

model = solph.Model(es)

model.solve("gurobi")

results = model.results()
17.959 <module>  profile_solph.py:1
├─ 6.740 Model.__init__  oemof\solph\_models.py:371
│  └─ 6.740 Model.__init__  oemof\solph\_models.py:91
│     └─ 6.728 Model._construct  oemof\solph\_models.py:164
│        ├─ 3.035 Model._add_child_blocks  oemof\solph\_models.py:190
│        │  ├─ 1.782 BusBlock._create  oemof\solph\buses\_bus.py:108
│        │  └─ 1.247 ConverterBlock._create  oemof\solph\components\_converter.py:220
│        │     ├─ 0.825 ConverterBlock.__setattr__  pyomo\core\base\block.py:551
│        │     │  └─ 0.825 ConverterBlock.add_component  pyomo\core\base\block.py:935
│        │     │     └─ 0.825 BuildAction.construct  pyomo\core\base\action.py:60
│        │     │        └─ 0.825 _input_output_relation  oemof\solph\components\_converter.py:255     
│        │     │           └─ 0.577 IndexedConstraint.add  pyomo\core\base\constraint.py:1019
│        │     └─ 0.305 IndexedConstraint.__init__  pyomo\core\base\constraint.py:742
│        ├─ 2.217 Model._add_objective  oemof\solph\_models.py:205
│        │  ├─ 1.473 SimpleFlowBlock._objective_expression  oemof\solph\flows\_simple_flow_block.py:407
│        │  │  ├─ 0.603 Series.__getitem__  pandas\core\series.py:1086
│        │  │  ├─ 0.362 [self]  oemof\solph\flows\_simple_flow_block.py
│        │  │  └─ 0.302 IndexedVar.__getitem__  pyomo\core\base\var.py:1049
│        │  └─ 0.695 Model.component_data_objects  pyomo\core\base\block.py:1543
│        │     └─ 0.671 Model._component_data_itervalues  pyomo\core\base\block.py:1451
│        │        ├─ 0.255 <genexpr>  pyomo\core\base\set.py:4142
│        │        └─ 0.248 IndexedVar.__getitem__  pyomo\core\base\var.py:1049
│        ├─ 1.145 Model._add_parent_block_variables  oemof\solph\_models.py:490
│        │  ├─ 0.515 Model.__setattr__  pyomo\core\base\block.py:551
│        │  └─ 0.243 IndexedVar.__getitem__  pyomo\core\base\var.py:1049
│        └─ 0.332 Model._add_parent_block_sets  oemof\solph\_models.py:401
│           └─ 0.318 Model.__setattr__  pyomo\core\base\block.py:551
├─ 4.445 Model.solve  oemof\solph\_models.py:238
│  └─ 4.442 GUROBISHELL.solve  pyomo\opt\base\solvers.py:531
│     ├─ 2.279 GUROBISHELL._presolve  pyomo\solvers\plugins\solvers\GUROBI.py:214
│     │  └─ 2.279 GUROBISHELL._presolve  pyomo\opt\solver\shellcmd.py:213
│     │     └─ 2.278 GUROBISHELL._presolve  pyomo\opt\base\solvers.py:688
│     │        └─ 2.278 GUROBISHELL._convert_problem  pyomo\opt\base\solvers.py:755
│     │           └─ 2.278 convert_problem  pyomo\opt\base\convert.py:25
│     │              └─ 2.262 PyomoMIPConverter.apply  pyomo\solvers\plugins\converter\model.py:44    
│     │                 └─ 2.258 Model.write  pyomo\core\base\block.py:1883
│     │                    └─ 2.257 LPWriter.__call__  pyomo\repn\plugins\lp_writer.py:193
│     │                       └─ 2.250 LPWriter.write  pyomo\repn\plugins\lp_writer.py:211
│     │                          └─ 2.250 _LPWriter_impl.write  pyomo\repn\plugins\lp_writer.py:250   
│     │                             ├─ 0.830 QuadraticRepnVisitor.walk_expression  pyomo\core\expr\visitor.py:257
│     │                             ├─ 0.466 _LPWriter_impl.write_expression  pyomo\repn\plugins\lp_writer.py:576
│     │                             ├─ 0.270 [self]  pyomo\repn\plugins\lp_writer.py
│     │                             └─ 0.204 Model.component_data_objects  pyomo\core\base\block.py:1543
│     ├─ 1.153 GUROBISHELL._apply_solver  pyomo\opt\solver\shellcmd.py:245
│     │  └─ 1.146 GUROBISHELL._execute_command  pyomo\opt\solver\shellcmd.py:308
│     │     ├─ 0.953 run  subprocess.py:506
│     │     └─ 0.191 TeeStream.__exit__  pyomo\common\tee.py:429
│     ├─ 0.759 ModelSolutions.load_from  pyomo\core\base\PyomoModel.py:186
│     └─ 0.203 GUROBISHELL._postsolve  pyomo\solvers\plugins\solvers\GUROBI.py:562
├─ 3.785 <module>  oemof\solph\__init__.py:1
│  ├─ 1.758 <module>  oemof\solph\buses\__init__.py:1
│  ├─ 1.040 <module>  oemof\solph\components\__init__.py:1
│  │  └─ 1.025 <module>  oemof\solph\components\experimental\__init__.py:1
│  │     ├─ 0.704 <module>  oemof\solph\components\experimental\_generic_caes.py:1
│  │     └─ 0.314 <module>  oemof\solph\components\experimental\_sink_dsm.py:1
│  └─ 0.943 <module>  oemof\solph\processing.py:1
│     └─ 0.940 <module>  pandas\__init__.py:1
│        └─ 0.620 <module>  pandas\core\api.py:1
└─ 2.968 Model.results  oemof\solph\_models.py:232
   └─ 2.968 results  oemof\solph\processing.py:188
      ├─ 2.684 create_dataframe  oemof\solph\processing.py:72
      │  ├─ 0.778 Series.map  pandas\core\series.py:4611
      │  ├─ 0.649 IndexedVar.__str__  pyomo\core\base\component.py:612
      │  ├─ 0.357 <listcomp>  oemof\solph\processing.py:81
      │  └─ 0.209 DataFrame.sort_values  pandas\core\frame.py:6984
      └─ 0.227 <dictcomp>  oemof\solph\processing.py:235

linopy

import pandas as pd
import linopy


m = linopy.Model()

hours = pd.Index(range(8760))

flow_names = ["source_b1", "b1_converter", "converter_b2", "b2_sink"]

flows = m.add_variables(lower=0, name="all-flows", coords=[hours, flow_names])

m.add_constraints(flows.loc[:, "b2_sink"] == 1, name="fixed demand")
m.add_constraints(flows.loc[:, "b1_converter"] == flows.loc[:, "converter_b2"], name='converter_efficiency')

m.add_constraints(flows.loc[:, "source_b1"] == flows.loc[: "b1_converter"], name="b1 balance")
m.add_constraints(flows.loc[:, "converter_b2"] == flows.loc[:, "b2_sink"], name="b2 balance")

variable_costs = 1
m.objective = variable_costs * flows.loc[:, "source_b1"]

m.solve(solver_name="highs")
m.objective.value

pd.pivot_table(index="dim_0", columns="dim_1", data=m.solution.to_dataframe())
# reads all variables into a dataframe
8.739 <module>  profile_linopy.py:1
├─ 3.901 Model.solve  linopy\model.py:928
│  └─ 3.781 run_gurobi  linopy\solvers.py:568
│     ├─ 1.603 to_file  linopy\io.py:261
│     │  ├─ 1.200 constraints_to_file  linopy\io.py:109
│     │  │  ├─ 0.551 [self]  linopy\io.py
│     │  │  ├─ 0.361 handle_batch  linopy\io.py:31
│     │  │  │  ├─ 0.201 [self]  linopy\io.py
│     │  │  │  └─ 0.145 IncrementalEncoder.encode  encodings\cp1252.py:18
│     │  │  └─ 0.202 Constraint.flat  linopy\constraints.py:456
│     │  └─ 0.320 bounds_to_file  linopy\io.py:170
│     ├─ 1.267 safe_get_solution  linopy\solvers.py:111
│     └─ 0.910 [self]  linopy\solvers.py
├─ 2.572 <module>  linopy\__init__.py:1
│  ├─ 1.989 <module>  linopy\monkey_patch_xarray.py:1
│  │  ├─ 1.584 <module>  linopy\expressions.py:1
│  │  │  ├─ 1.409 <module>  scipy\sparse\__init__.py:1
│  │  │  │  ├─ 1.130 <module>  scipy\sparse\csgraph\__init__.py:1
│  │  │  │  │  ├─ 0.839 <module>  scipy\sparse\csgraph\_laplacian.py:1
│  │  │  │  │  └─ 0.247 [self]  scipy\sparse\csgraph\__init__.py
│  │  │  │  └─ 0.098 <module>  scipy\sparse\_base.py:1
│  │  │  │     └─ 0.090 <module>  scipy\_lib\_util.py:1
│  │  │  │        └─ 0.088 <module>  scipy\_lib\_array_api.py:1
│  │  │  └─ 0.121 <module>  scipy\__init__.py:1
│  │  └─ 0.399 <module>  xarray\__init__.py:1
│  └─ 0.568 <module>  linopy\model.py:1
├─ 1.242 <module>  pandas\__init__.py:1
│  ├─ 0.711 <module>  pandas\core\api.py:1
│  └─ 0.292 <module>  numpy\__init__.py:1
├─ 0.624 Model.__init__  linopy\model.py:83
└─ 0.263 Variable.__eq__  linopy\variables.py:348
@fwitte
Copy link
Member Author

fwitte commented Mar 27, 2024

Fast model building would also be beneficial for rolling horizon dispatch planning! I will try to make a minimal implementation into solph and open a PR.

@fwitte
Copy link
Member Author

fwitte commented Apr 9, 2024

Here is a first implementation draft: https://github.com/oemof/oemof-solph/blob/feature/add-linopy-backend/solph-with-linopy-concept.py

It is not beautiful, but it may be a good point to start a discussion on if and how to move forward.

@metab0t
Copy link

metab0t commented May 27, 2024

Hello, @fwitte!

I am the author of PyOptInterface, an efficient modeling interface for mathematical optimization in Python. Its performance is quite competitive compared with existing solutions (faster than vendored Python bindings of some optimizers).

PyOptInterface provides user-friendly expression-based API to formulate problems with different optimizers, and the lightweight handles of variables and constraints can be stored freely in built-in multidimensional container or Numpy ndarray.

By the way, PyOptInterface supports efficient incremental modification of the model (adding/removing variables or constraints, modifying the coefficient and RHS of linear constraints, change the objective, etc). You can refer to the documentation of deleting variables, deleting constraints, modifying constraints and modifying objective. This might be very useful for rolling dispatch problems.

As far as I know, the time to modify and solve in linopy is quite expensive because it needs to re-export the model again.

If you are concerned with the performance to construct models, welcome to take a look at PyOptInterface and try it!

I have prepared a benchmark between linopy and PyOptInterface where we construct a model, then modify and solve it in iterations.

import time

import pandas as pd
import linopy

import pyoptinterface as poi
from pyoptinterface import gurobi


def linopy_main(N, K):
    m = linopy.Model()

    hours = pd.Index(range(N))

    flow_names = ["source_b1", "b1_converter", "converter_b2", "b2_sink"]

    flows = m.add_variables(lower=0, name="all-flows", coords=[hours, flow_names])

    fixed_demand_con = m.add_constraints(
        flows.loc[:, "b2_sink"] == 1, name="fixed demand"
    )
    m.add_constraints(
        flows.loc[:, "b1_converter"] == flows.loc[:, "converter_b2"],
        name="converter_efficiency",
    )

    m.add_constraints(
        flows.loc[:, "source_b1"] == flows.loc[:"b1_converter"], name="b1 balance"
    )
    m.add_constraints(
        flows.loc[:, "converter_b2"] == flows.loc[:, "b2_sink"], name="b2 balance"
    )

    variable_costs = 1
    m.objective = variable_costs * flows.loc[:, "source_b1"]

    m.solve(solver_name="gurobi")

    for k in range(K):
        fixed_demand_con.rhs.loc[k] = 2.0
        m.solve(solver_name="gurobi")


def poi_main(N, K):
    m = gurobi.Model()

    hours = range(N)

    flow_names = ["source_b1", "b1_converter", "converter_b2", "b2_sink"]

    flows = m.add_variables(hours, flow_names, lb=0, name="all-flows")

    def fixed_demand(h):
        return m.add_linear_constraint(
            flows[h, "b2_sink"], poi.Eq, 1.0, name=f"fixed demand {h}"
        )

    fixed_demand_con = poi.make_tupledict(hours, rule=fixed_demand)

    def converter_efficiency(h):
        return m.add_linear_constraint(
            flows[h, "b1_converter"] - flows[h, "converter_b2"],
            poi.Eq,
            0.0,
            name=f"converter efficiency {h}",
        )

    converter_efficiency_con = poi.make_tupledict(hours, rule=converter_efficiency)

    def b1_balance(h):
        return m.add_linear_constraint(
            flows[h, "source_b1"] - flows[h, "b1_converter"],
            poi.Eq,
            0.0,
            name=f"b1 balance {h}",
        )

    b1_balance_con = poi.make_tupledict(hours, rule=b1_balance)

    def b2_balance(h):
        return m.add_linear_constraint(
            flows[h, "converter_b2"] - flows[h, "b2_sink"],
            poi.Eq,
            0.0,
            name=f"b2 balance {h}",
        )

    b2_balance_con = poi.make_tupledict(hours, rule=b2_balance)

    variable_costs = 1
    objective = poi.quicksum(variable_costs * flows[h, "source_b1"] for h in hours)

    m.set_objective(objective)

    m.optimize()

    for k in range(K):
        m.set_normalized_rhs(fixed_demand_con[k], 2.0)
        m.optimize()


def main(N, K):
    linopy_times = {}
    for k in range(K):
        start = time.time()
        linopy_main(N, k)
        end = time.time()
        linopy_times[k] = end - start

    poi_times = {}
    for k in range(K):
        start = time.time()
        poi_main(N, k)
        end = time.time()
        poi_times[k] = end - start

    for k in range(K):
        print(f"Run ({N}, {k}):")
        linopy_time = linopy_times[k]
        poi_time = poi_times[k]
        print(f"\tlinopy: {linopy_time:.2f} seconds")
        print(f"\tpoi: {poi_time:.2f} seconds")


if __name__ == "__main__":
    main(8760, 10)

On my computer, the output is:

Run (8760, 0):
        linopy: 0.75 seconds
        poi: 0.20 seconds
Run (8760, 1):
        linopy: 1.31 seconds
        poi: 0.16 seconds
Run (8760, 2):
        linopy: 1.94 seconds
        poi: 0.15 seconds
Run (8760, 3):
        linopy: 2.58 seconds
        poi: 0.16 seconds
Run (8760, 4):
        linopy: 3.19 seconds
        poi: 0.18 seconds
Run (8760, 5):
        linopy: 3.89 seconds
        poi: 0.18 seconds
Run (8760, 6):
        linopy: 4.49 seconds
        poi: 0.19 seconds
Run (8760, 7):
        linopy: 5.18 seconds
        poi: 0.20 seconds
Run (8760, 8):
        linopy: 5.90 seconds
        poi: 0.21 seconds
Run (8760, 9):
        linopy: 6.85 seconds
        poi: 0.26 seconds

As the result shows, the one-shot time of PyOptInterface is nearly 1/4-1/3 of linopy. The time of linopy increases linearly as the number of "modify and solve" iterations K grows, while PyOptInterface remains efficient because it modifies the changed part in-place without rebuilding the entire model.

@p-snft
Copy link
Member

p-snft commented May 29, 2024

I checked other boundary conditions:

  • The license seems to be fine. The package states that it is MPL 2.0, which is a weak copyleft license. Thus, linking or importing has no legal implications for our package per se. Only distributions that include PyOptInterface also have to include the PyOptInterface code and the MPL.
  • At least on Linux machines (esp. for containers as used in automated pipelines), HiGHS is rather easy to install. I'd even consider it an advance compared to Pyomo with cbc.
  • I did not check code quality, yet. (PyOptInterface being a relatively new one-person project isn't really convincing.)

@metab0t
Copy link

metab0t commented May 29, 2024

Thanks for your feedback!

For the third concern, the design of PyOptInterface is a minimalistic and thin wrapper of C API with no magic, so its implementations are essentially calling corresponding C API. You can look at https://github.com/metab0t/PyOptInterface/blob/master/lib/gurobi_model.cpp as an example. The logic of each function is quite clear and explicit.

On the Python side, we only provide one built-in container type tupledict following the similar semantics of tupledict in gurobipy. It is up to the users to use high-level structures to store these variables and constraints freely, such as numpy.ndarray or pandas.DataFrame. There is no layered abstractions to make things more complex. (The tupledict we provide is just a Python dict with extra methods)

As for the maturity of PyOptInterface, some academic researchers I know (including myself) have been using it in research and engineering on power systems and have positive experience. The minimalistic design of PyOptInterface also significantly reduces its maintenance burden because the C API of optimizers are very stable even across major versions.

Of course, I am looking forward to a bigger community of users and contributors as it gets more adoptions.

@p-snft
Copy link
Member

p-snft commented May 29, 2024

The third point was no objection, rather a reason why checking code quality might be very important. Indeed, even on the long run, a slim one-person-project can be better maintained than a big, established one.

Running the benchmark you posted on my local machine with HiGHS as a solver resulted in the following output:

Run (8760, 0):
	linopy: 1.05 seconds
	poi: 285.46 seconds
Run (8760, 1):
	linopy: 1.59 seconds
	poi: 266.68 seconds
Run (8760, 2):
	linopy: 2.74 seconds
	poi: 263.65 seconds
Run (8760, 3):
	linopy: 4.16 seconds
	poi: 219.76 seconds
Run (8760, 4):
	linopy: 5.38 seconds
	poi: 203.12 seconds
Run (8760, 5):
	linopy: 6.73 seconds
	poi: 202.25 seconds
Run (8760, 6):
	linopy: 8.00 seconds
	poi: 139.61 seconds
Run (8760, 7):
	linopy: 9.15 seconds
	poi: 136.08 seconds
Run (8760, 8):
	linopy: 10.21 seconds
	poi: 138.26 seconds
Run (8760, 9):
	linopy: 11.47 seconds
	poi: 140.75 seconds

At least without adjusting the setup, linopy outperforms poi by far. Interestingly, the time to solve the poi model decreases for higher number of iterations. It feels like something strange is happening there.

@p-snft p-snft changed the title Implement linopy as alternative back-end to pyomo Implement alternative back-end to Pyomo (e.g. linopy) May 29, 2024
@metab0t
Copy link

metab0t commented May 29, 2024

I can reproduce this result. Thank you very much for the test on HiGHS.

I have figured out that there are four HiGHS API functions that slow down PyOptInterface significantly. I have a local fixed version that makes HiGHS fast again, and I will release a new version soon.

@metab0t
Copy link

metab0t commented May 29, 2024

@p-snft I have pushed a new release v0.2.2. Use pip install -U pyoptinterface to test it and it is fast again!

This is the result of HiGHS on my computer:

Run (8760, 0):
        linopy: 0.76 seconds
        poi: 0.32 seconds
Run (8760, 1):
        linopy: 0.99 seconds
        poi: 0.24 seconds
Run (8760, 2):
        linopy: 1.27 seconds
        poi: 0.27 seconds
Run (8760, 3):
        linopy: 1.56 seconds
        poi: 0.24 seconds
Run (8760, 4):
        linopy: 1.84 seconds
        poi: 0.25 seconds
Run (8760, 5):
        linopy: 2.12 seconds
        poi: 0.27 seconds
Run (8760, 6):
        linopy: 2.48 seconds
        poi: 0.28 seconds
Run (8760, 7):
        linopy: 2.64 seconds
        poi: 0.29 seconds
Run (8760, 8):
        linopy: 3.03 seconds
        poi: 0.30 seconds
Run (8760, 9):
        linopy: 3.24 seconds
        poi: 0.31 seconds

The first time is a little slower because PyOptInterface needs to load the dynamic library of HiGHS in the initialization step.

@metab0t
Copy link

metab0t commented May 30, 2024

The $O(N^2)$ performance issue of HiGHS is fixed in ERGO-Code/HiGHS#1782 by the upstream.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants