Ask performance tips

[Zhanwei-Liu]

Dear developers of PyOptInterface,

I recently used PyOptInterface to replace Pyomo for my energy expansion model called PREP-SHOT. However, I have not observed a significant improvement in compiling speed. For a model with over 6 million variables, the compiling time decreased from 199.89 seconds with Pyomo to 144.08 seconds with PyOptInterface on a Mac OS Apple M1 Max with 10 cores. Do you have any suggestions for further improving the compiling speed? I suspect I might not be using the most efficient coding practices.

Thank you very much.

Best regards,

[metab0t]

Thanks for trying out our project.
Do you have any benchmark script that I can run directly to test the performance? I see the run.py in the root directory but am not sure how to use it.

[Zhanwei-Liu]

Thank you for your reply!

Please follow the procedure below to download the benchmark script. I apologize for any inconvenience, but I want to use an example that approximates a real-world and challenging scenario.

git clone -b PyOptInterface https://github.com/PREP-NexT/PREP-SHOT.git
cd PREP-SHOT
conda env create -f prep-shot.yml
conda activate prep-shot
python run.py

You can modify the solving configuration in the solve_parameters section of config.json. All mathematical modeling is included in prepshot/model.py and prepshot/rules.py.

Thank you very much!

[metab0t]

I fix some performance issues in this commit.

The main issue is to avoid using code like this:

expr = poi.quicksum(...)
lhs = x - expr

Use in-place operation instead:

expr = poi.quicksum(...)
lhs = expr
lhs -= x

Because the +-*/ operators of ExprBuilder will copy its content, please use +=,-=,*=,/= these in-place operators on ExprBuilder to construct the big expression efficiently.

On my computer, the model building time decreases from 374s to 248s (with Gurobi optimizer). You can try it and share your results.

[Zhanwei-Liu]

Thank you for your insights and commit!

I'm not sure why, but when I used the merged codes to re-run the model, the model building time only slightly decreased (using the Gurobi optimizer). I tried it twice: the time decreased from 143.45 seconds to 143.18 seconds in the first run and from 141.94 seconds to 141.07 seconds in the second run.

Thank you!

[metab0t]

I think the problem is that your measured "problem compiling" time is only a small part of the time spent by pyomo. The solver.solve() in Pyomo is the most time-consuming part where pyomo passes its internal model to the solver. Meanwhile, model.optimize() is very fast in PyOptInterface because the model has been constructed in the optimizer.

You should set the time limit to 0.0 seconds and measure the total time.

[Zhanwei-Liu]

Thank you for your suggestion.

You are correct! I set the time limit to 0.0 seconds and measured the total time twice. Pyomo took 474.3 and 478.0 seconds, while PyOptInterface took 142.2 and 141.8 seconds.

For your reference, I also rewrote the model-building part using JuMP based on PyJulia. JuMP took 71.3 and 69.9 seconds. However, when I switched from PyJulia to PythonJulia, the time increased to 102.0 and 105.2 seconds. To conclude, I should select PyOptInterface as a dependency for my model, as it is more appropriate from a long-term perspective.

I believe my question is resolved. Thank you very much.

[metab0t]

Your speed is so fast!
The indexing speed of tupledict is much slower than containers in JuMP due to the internal weakness of Python. tupledict is more similar to the SparseAxisArray in JuMP.
It would be great to implement DenseAxisArray in Python with numpy, but it is very difficult to be faster than Julia in this aspect.
Besides, I assume that the memory consumption of POI is also less than Pyomo?

[metab0t]

@Zhanwei-Liu Would you like to test my latest commit febd996 and change your local tupledict.py to test its influence on the performance?
I think it may make make_tupledict faster.

[Zhanwei-Liu]

Thank you very much!

I tested the performance after updating POI. The speed improvement is minimal, with time decreasing from 141 seconds to 138.62 and 137.00 seconds in two experiments. I suspect the generation of expressions is the main bottleneck.

Additionally, the memory test results were unexpected. I am not very familiar with how to test memory usage, so I used the following script for the test:

import psutil

def timer(func):
    """A decorator that times a function.

    Parameters
    ----------
    func : function
        The function to be timed.

    Returns
    -------
    function
        The decorated function.
    """
    @functools.wraps(func)
    def wrapper(*args, **kwargs):
        process = psutil.Process(os.getpid())
        start_memory = process.memory_info().rss
        start_time = time.perf_counter()
        value = func(*args, **kwargs)
        end_memory = process.memory_info().rss
        end_time = time.perf_counter()
        run_time = end_time - start_time
        memory_used = (end_memory - start_memory) / 1024 / 1024  # MB
        logging.info(
            "Finished %s in %.2f seds", repr(func.__name__),
            run_time
        )
        logging.info(
            "Memory used %s in %.2f MB", repr(func.__name__),
            memory_used
        )

        return value

    return wrapper

The test results in two experiments:

POI: 6065.56 MB and 6084.22 MB
Pyomo: 3746.44 MB and 3742.89 MB

[metab0t]

We should measure the peak memory usage in the solution process. Can you try the following code:

stop_monitoring = False
# List to store memory usage values
memory_values = []
if __name__ == "__main__":
    # Function to monitor memory usage
    def monitor_memory_usage(interval=0.1):
        global stop_monitoring
        global memory_values
        process = psutil.Process()
        while not stop_monitoring:
            mem_info = process.memory_info()
            memory_values.append(mem_info.rss / 1024**2)  # Store memory in MB
            time.sleep(interval)

    # Start monitoring memory usage in a separate thread
    monitor_thread = threading.Thread(target=monitor_memory_usage)
    monitor_thread.daemon = True  # Daemonize thread
    monitor_thread.start()

    main()

    # stop monitoring
    stop_monitoring = True
    monitor_thread.join()

    max_memory = max(memory_values)
    logging.info(f"Maximum memory usage: {max_memory:.2f} MB")

BTW, I do not know how to test Pyomo with the southeast_asia input data.

[Zhanwei-Liu]

Thank you!

I modified my timer decorator function by using a Python built-in library called tracemalloc to test peak memory usage due to its convenience. Below is the code I used:

import tracemalloc
def timer(func):
    """
    Decorator to log the start and end of a function, and how long it took to run.

    Args:
        func (function): The function to be decorated.

    Returns:
        function: The decorated function.
    """
    def wrapper(*args, **kwargs):
        # logging.info("Start solving model ...")
        start_time = time.time()
        tracemalloc.start()
        result = func(*args, **kwargs)
        _, peak = tracemalloc.get_traced_memory()
        tracemalloc.stop()
        run_time = time.time() - start_time
        logging.info(
            "Finished %s in %.2f seds", repr(func.__name__),
            run_time
        )
        logging.info(
            "Memory used %s in %.2f MB", repr(func.__name__),
            peak / 10**6
        )
        return result
    return wrapper

The test results are as follows:

• POI: 2004.31 MB
• Pyomo: 16882.72 MB

Clearly,Pyomo uses over 8 times the memory of POI for this example. It's exciting!

For your reference, I have uploaded a new branch called Pyomo. You can reproduce my results by cloning the Pyomo and PyOptInterface branches if you want.

[metab0t]

The improvement on memory usage is quite significant! Thanks for your work and detailed benchmark.