Linear Optimization in time-critical Control Loop

13 May, 2023

Problem

Given a start and endpoint, find the eager and lazy path from start to end.

Model

The visualization translates to a simple rhombus (German: Raute) where the enclosed area is to be maximized.

The task is to find the corner points A, B, C, D, of the Raute.

There is a twist, however. The points must be contained within a “valid” area, whereas start and endpoint may lie outside.

Without this twist, the problem would be trivial to solve with a simple linear equation system.

The problem still can be solved easily with simple linear math, but due to the twist it would require conditional calculations, depending whether the start/endpoint lie within or outside the valid area.

With a LP (linear programming) model a closed formulation without conditional branching is possible:

Parameters

$$ S=(t_S, Y_S)\\ E=(t_E, Y_E)\\ P_0 > 0\\ C_1 = Y_S - P_0 t_S \\ C_2 = Y_E - P_0 t_E \\ $$

Objective

$$ \max \ t_C - t_A \\ $$

Variables

$$ A=(t_A, Y_A)\\ B=(t_B, Y_B)\\ C=(t_C, Y_C)\\ D=(t_D, Y_D)\\ $$

Feasability

$$ C_1 \ge C_2\\ $$

Constraints

Allowed Area

$$ Y_{min} \le P_0 t_A + C_1 \le Y_{max}\\ Y_{min} \le P_0 t_C + C_1 \le Y_{max}\\ Y_{min} \le P_0 t_B + C_2 \le Y_{max}\\ Y_{min} \le P_0 t_D + C_2 \le Y_{max}\\ $$

Start/End

$$ Y_{S} \le P_0 t_A + C_1 \le Y_{E}\\ Y_{S} \le P_0 t_C + C_1 \le Y_{E}\\ Y_{S} \le P_0 t_B + C_2 \le Y_{E}\\ Y_{S} \le P_0 t_D + C_2 \le Y_{E}\\ $$

Raute

$$ Y_A = Y_B \\ Y_C = Y_D \\ Y_A = P_0 t_A + C_1 \\ Y_C = P_0 t_C + C_1 \\ Y_B = P_0 t_B + C_2 \\ Y_D = P_0 t_D + C_2 \\ $$

With Pyomo it is easy to translate the LP model into Python code. The model solves in a split second.

Result

Start/End outside valid area

Solution 35.0
Number of constraints : 14
Number of variables : 8

{'Y': {'A': 50.0, 'B': 50.0, 'C': 400.0, 'D': 400.0},
 't': {'A': 14.0, 'B': 55.0, 'C': 49.0, 'D': 90.0}}

Start/End within valid area

Solution 28.0
Number of constraints : 14
Number of variables : 8

{'Y': {'A': 70.0, 'B': 70.0, 'C': 350.0, 'D': 350.0},
 't': {'A': 10.0, 'B': 72.0, 'C': 38.0, 'D': 100.0}}

High Performance Solution

Pyomo is very capable and flexible for any modelling needs, but performance-wise we can do better.

Before the solver can start working on the problem, Pyomo needs to translate the model into an LP file format which will then be read by the (CBC) solver. This indirection is an “ok” worflow for quick iteration and model building.

However, if for example the optimization needs to run in a critical real-time control loop this might be inacceptable.

A better way is to directly talk to the solver API. Here Rust and good_lp provide an interesting option.

Translating the model into Rust allows to directly bind the shared libraries of CBC and avoid any intermediate translation steps.

Rust Implementation

Translating the model into Rust is straight forward for such a simple problem:

    // Objective
    let solution = vars.maximise(t_C - t_A)
        .using(good_lp::default_solver)
        .with(constraint!(Ymin <= E_A ))
        .with(constraint!(Ymin <= E_B ))
        .with(constraint!(Ymin <= E_C ))
        .with(constraint!(Ymin <= E_D ))
        .with(constraint!(E_A <= Ymax))
        .with(constraint!(E_B <= Ymax))
        .with(constraint!(E_C <= Ymax))
        .with(constraint!(E_D <= Ymax))

        .with(constraint!(S.Y <= E_A ))
        .with(constraint!(S.Y <= E_B ))
        .with(constraint!(S.Y <= E_C ))
        .with(constraint!(S.Y <= E_D ))
        .with(constraint!(E_A <= E.Y))
        .with(constraint!(E_B <= E.Y))
        .with(constraint!(E_C <= E.Y))
        .with(constraint!(E_D <= E.Y))

        .with(constraint!(E_A == E_B))
        .with(constraint!(E_C == E_D))

        .with(constraint!(E_A == P * t_A + Ce))
        .with(constraint!(E_C == P * t_C + Ce))
        .with(constraint!(E_B == P * t_B + Cl))
        .with(constraint!(E_D == P * t_D + Cl))
        .solve()?;

Comparison

The “rustified” model is more than 50 times faster then the Pyomo version! This is impressive, because the actual solver code is in both cases highly optimized CBC C++ code.

Solving the model with Pyomo is not the bottleneck, CBC is fast once it has read the LP file. But generating the LP file and reading it into CBC adds a lot of overhead.

Pyomo

Rust

Conclusion

To model a LP/MIP problem Pyomo is a very attractive solution. In any case it is much more powerful and flexible than Rust good_lp.

Combining Pyomo’s strengths with the raw performance of Rust and direct library binding results in a very competitive and capable optimization toolchain for LP and MIP problems.

#python #optimization #pyomo #rust